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guidance scheme with appropriate time lags is employed to model the essential 
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bances. The results indicate that for certain regions of parameters, limit cycle 
oscillations may develop which could compromise system stability and safety 
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I. INTRODUCTION 
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A. PROBLEM STATEMENT 

The problem of motion stability of ships and other marine vehicles has been 
the subject of extensive studies in the past (Comstock, 1977). Most of these 
studies are with regards the ship’s directional stability in open waters under 
open loop conditions. It is well known that in restricted waters such as canals 
or rivers, although it is possible to have positional stability in principle, in 
reality it is never the case. This is due to the destabilizing effects of the bank 
suction forces and moments. These act always in a destabilizing fashion; i.e., 
after a small initial disturbance they produce a force or a moment which tends 
to increase the ship’s path deviation from nominal (Comstock, 1977). Open 
loop conditions, important as they are, rarely represent real life applications. 
Ships traveling in canals are under closed loop control, typically provided 
through the helmsman or an autopilot. It becomes then important to assess 
the stability of the system under finite disturbances and incorporating some 
modeling of closed loop operations. 



B. OBJECTIVES AND OUTLINE 

This thesis utilizes both linear and nonlinear techniques in order to analyze 
the problem of closed loop dynamic stability of ships in canals. Ship modeling 
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is discussed in Chapter II and it follows standard ship maneuvering equations 
(Clayton and Bishop, 1982) augmented by a model for bank suction effects 
(Beck, 1976). A Mariner class ship is selected as a model because of the avail- 
ability of data for its hydrodynamic coefficients (Comstock, 1977), (Bernitsas 
and Kekridis, 1984). A heading control law based on Nomoto’s first-order 
model, coupled with a line of sight guidance law is chosen in order to model 
the fundamental behavior of closed loop conditions. Since, in practice control 
actions are hardly ever instantaneous, appropriate time lags are introduced 
in the feedback. These are modeled using the techniques outlined in (Venne, 
1992). Linear eigenvalue analysis is performed in Chapter III in order to 
reveal parametric stability boundaries. It is established that loss of stability 
occurs when a pair of complex conjugate eigenvalues crosses the imaginary axis 
which results in periodic solutions or limit cycles (Guckenheimer and Holmes, 
1983). A nonlinear analysis based on Hopf bifurcation methods (Chow and 
Mallet-Paret, 1977) and (Hassard and Wan, 1978) is performed in Chapter 
IV. Conclusions derived from this study and some recommendations for fur- 
ther extensions are discussed in Chapter V. The basic programs that were used 
to produce the results are included in the Appendix. 
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II. PROBLEM FORMULATION 



A. SHIP DYNAMICS 

Restricting our attention to the horizontal plane, the mathematical model 
consists of the nonlinear sway (translational motion parallel to the vehicle’s 
longitudinal axis) and yaw (rotational motion about the vertical axis ) equa- 
tions of motion. If we consider a moving Cartesian coordinate frame located 
at the geometric center of the vehicle, Newton’s equations of motion are: 



m{v + ur + XGT) = Y, 


(1) 


IzT + mxciv -1- ur) = N , 


(2) 



where v and r are the relative sway and yaw velocities of the moving vehicle 
with respect to the water, m is the vehicle’s mass, Iz is its moment of in- 
ertia with respect to the vertical axis, and u is the constant forward speed. 
Since all quantities will be considered as dimensionless in this study, in the 
following we consider u to be equal to one. xq is the longitudinal position 
of the ship’s center of gravity with respect to its centroid, and Y,N represent 
the total excitation sway force and yaw moment respectively. Following stan- 
dard ship maneuvering assumptions, these forces can be expressed as the sum 
of quadratic drag terms and first-order memoryless polynomials in v and r, 
which represent added mass and damping due to water. In this way they can 
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be expressed by: 



N 



Y 




( 3 ) 



(4) 



where Ya, Na represent the partial derivatives with respect to the indicated 
variable a and 6 is the effective rudder angle. Y{ip,y), N{'ip,y) represent the 
interaction/proximity forces and moments that arise due to the presence of 
the canal, and shallow water effects. 

The resulting nonlinear differential equations can be non-dimensionalized 
with respect to the constant forward speed of the ship, u and its length, 1. 
The dimensionless time variable is then equal to tu/l. A standard inertial sys- 
tem (x,y) is introduced here where the x-axis points in the assumed nominal 
straight line path and the y-axis is the distance from this nominal path, see 
Figure 1. We assume that the nominal straight line path is along the center- 
line of the canal. In cases where the concept of a geometric centerline is not 
applicable, we assume that the nominal path is along the zero bank suction 
location. This is consistent with recommended navigation practices that are 
currently in use. 

The complete model is presented by the following equations: 
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Canal centerline <■ 



Xd 



■> 




Figure 1: Vehicle geometry and definitions of symbols 

+ {Yr - m)r + F (^, y) + Ys6 , (5) 

— {Ny - mxa)v + {h — Nr)r = NyV + -NyyyV^ + -NyrrVr‘^+ , 

6 2 

-Nryyrv^ + {Nr - mxa)r + Ni^p, y) + Ns6 , (6) 

ijj=r, (7) 

y = sim/; + V cos V' , (8) 



where the expressions for the ship’s yaw rate as well as its inertial position rate 
have been added, and xp is the local heading angle. The inter action/proximity 
forces and moments are expanded to include up to third order terms, 

Y {ip, y) = Y^xp + Yyy + ^Y^^^xp^ + Wyyyy^ , 

0 t) 

N {xp,y) = N-^xp + Nyy + —N-^^^xp^ 4- ^Nyyyy^ , 
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( 9 ) 

( 10 ) 



as explained in more detail in the folowing section. 



B. HYDRODYNAMIC COEFFICIENTS 

We chose a Mariner class ship as a representative model. Its hydrodynamic 
coefficients and geometric properties were taken from (Comstock, 1977) and 
(Bernitsas and Kekridis, 1984). Results from (Beck, 1976) were used in order 
to model the bank suction forces and moments. Typical results are shown in 
Figure 2. These show force and moment coefficients versus ship deviation (rj), 
and for canal width w = 0.4L. Force and moment coefficients are nondimen- 
sionalized with respect to the water density and the ship’s speed and length, 
as is customary in ship maneuvering. Since the suction force and moment 
must change their sign as either y oi x/j changes its sign, they must be modeled 
by odd power polynomials, as was done in the previous section. Numerical 
values for the coefficients Yy, Yyyy, Ny, and Nyyy can be found by curve fitting 
of the results of Figure 2. Using a depth to draft ratio of 1.9 we were able to 
estimate. 



Yy = 0 . 02 , 

^yyy ~ 0.468 , 

Ny = -0.0025, 

^yyy ^ • 

The value of was estimated by “discretizing” the ship in two segments. 
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FIGURE 13 

VARIATION OF SIDE FORCE AND MOMENT WITH 
WALE POSITION RATIO FOR THE MARINER 



w/L = .4 





Figure 2: Forces and moments due to canal [Beck (1976)] 



scale 

1.3 

1.9 
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at the bow and the stern, and calculating the suction forces on each part 
individually. Their resultant moment produced, 

= 0.01 . 

The value of was taken to be zero, due to lack of reliable data. The value 

of was estimated to be equal to —Yy = 0.014 which is true for motions along 

the centerline of the canal (Comstock, 1977). Again, due to lack of reliable 
data we took Y^-^^ = 0. 

C. CONTROL LAW 

The linearized set of equations (5) through (7) can be expressed in the 
following form: 

ij} = r , 

V = aiiip + ai2V + oi3r + a^y + bi6 , (11) 

r = a2iip + a22V + a2zr + a2AV + , 

where the coefficients Oy, bi are functions of the vehicle hydrodynamic coeffi- 
cients and geometric properties and are given below: 

an = — ^[(/z - Nr)Y.^ + {Yf - mxG)N^] , 

Uyr 

a \2 = ~jr ~ Nr)^v “ 1 “ (y^r j 

Uyr 

ai3 = - Ni){Yr - m) + (Y^ - mxG){Nr - mxG)] , 

Uyr 

014 = 7^[(^2 - ^r)^y + (^r ~ mXG)Ny] , 

Uyr 
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021 = j—liNy - mxG)y'^ + (m - Yy)N^] , 

U'ur 



vr 




vr 



h = -^[{Iz-Nr)Y6 + {Yr-mXG)Ns], 

Uyr 




vr 



Dyr = i'm- Yy){Iz - Nr) - {Ny - mxG)(Yr - mxG) • 

A control law that could model a human operator should not be based on 
feedback of the side slip velocity v. Instead, it is more likely that human 
operators will respond to changes in the ship’s heading angle -tp and rate of 
change of heading, r. Therefore, we choose to base our control law on Nomoto’s 
model, which follows by assuming negligible sway velocity v, 



f = ar cxjj b6 , 



( 12 ) 



where the coefficients are 



a 



023 , 



b = b2, 



c 



021 • 



A linear control law can be introduced in the form 



^0 = kii'ip - ‘ipc) + 



( 13 ) 
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where is the commanded heading angle and the gains k \ , k 2 are computed 
such that the closed loop system (7), (12), and (13) has the desired dynamics. 
The existence of the difference (■0 — 0c) in the control law (13) has the effect of 
trying to point the ship’s longitudinal axis towards the desired heading. The 
characteristic equation of the system is obtained from (7), (12), and (13) as 

— (a + hk 2 )s — (c + hk\) = 0 . 

If the desired characteristic equation is 

+ 2C,U)nS + = 0 , 



the control gains can be computed from 




k2 



o, 2 (^lo 
6 



n 



The natural frequency uJn and damping ratio C are selected based on general 
properties of second order systems. Relatively high values of and low values 
of C will correspond to a responsive operator, while the opposite is true for 
combinations of low and high 



Finally, in order to take into account the effect of rudder saturation, the 
commanded rudder angle is given by 

6 = 6sat tanh , (14) 

where (5o is the slope of 6 at the origin given by (13), and 6sat is the saturation 
limit on 5, typically around 0.4 radians. The hyperbolic tangent function is 
used instead of a hard saturation function, because of its differentiability. 
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D. GUIDANCE SCHEME 



Since the previous control law stabilizes the ship to any commanded head- 
ing angle, it must be coupled with an appropriate orientation guidance scheme 
to provide path keeping along the desired track. The simplest guidance scheme 
which models some fundamental aspects of helmsman behavior is a pure pur- 
suit guidance where the commanded heading angle equals the line of sight 
angle, 




as shown in Figure 1. The ship is located at (x,y) and attempts to point its 
longitudinal axis towards a target point which is located ahead of the ship 
on the reference path at a constant preview distance Xd- Pursuit guidance is 
achieved by commanding a heading angle ipc equal to the line of sight angle 
(15). 

The positional error information y in equation (15), is assumed to lag the 
actual error y by an amount of T seconds, in other words, 

y = y{t-T). (16) 

In this equation, the time lag T models the necessary time that it takes for 
the helmsman to process his path deviation and initiate appropriate corrective 
actions. 
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III. LINEAR ANALYSIS 



In this Chapter, we present a linearized analysis of the equations of motion. 
The purpose of this analysis is to assess stability or instability of the equations 
in response to small deviations from the straight line reference track. No 
attempt will be made here to characterize the transient response of the system. 
This can be accomplished through numerical simulations. 



A. COMBINED SYSTEM 

If we incorporate the interaction/proximity forces (9) and (10) in the ship 
dynamics model, equations (5) through (8), we end up with the combined 
system. 

Ip = r , 

{in T^)u {'^r mXQ^T — ~z'^wv'^ “b "b ”b 

1 6 1 2 2 

(Yj- — m)r + Y.^'tp + Yyy + —Y^p^ip^ + —Yyyyi/^ + Ys6 , 

® 1^1 

-{Ny - mXG)V + {Iz - Nr)r = NyV + -NyyyV^ + -Nyrrvr^+ (17) 

1 ® 1 ^ 

-NryyVv^ + {Nr - mxG)r + N^tp + NyX) + 

? 6 

'Z^yyyy^ ■*” > 

o 

y = sin ■ip +v cos . 

Study of the asymptotic properties of this system is the subject of this and 
the following chapters. 
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B. LINEARIZATION 



The linearized form of equations (17) is the following, 

Ip , 

(m — Yv)v — (y,; — mxo)r = YyV + {Yr — m)r+ 

+ Vyj/ + 

-{Ni, - mxG)v + {h- Nr)r = NyV -{■ {Nr - mxG)r+ ^ 

N-^ip + NyD + Ns6 , 
y = 'ip + V . 

The rudder angle 6 has one of the forms that are developed below. The time 
delay operator can be expressed in terms of its Taylor series expansion, 

y{t-T)=y-Ty + ]-T^y - + • • • . (19) 

z o 

Practical computations can be performed by truncating equation (19) to first, 
second, or third order. 

1. First order approximation in y 
We have, 

y{t -T) = y - Ty , 



or 

y{t - T) = y - T{-tp + v) . 



In its linear form, 



and, therefore. 





J 



6 = 6q. 



Furthermore, in its linear form equation (15) can be written as 

Xd 



( 20 ) 
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If we incorporate (20) into (13), we derive the linearized first order approxi- 



mation in 6, 



. , , fciT . . 

6 = ki-ij} + k 2 T H y (■0 + v) . 

Xd Xd 



( 21 ) 



2. Second order approximation in y 

Keeping the second order terms in y{t — T) we get, 

y{t-T) = y-Ty + ]^T‘^y . (22) 

If we incorporate (22) into (13) along with the linear equations 6 — 6q and 



= -y{t - T)/xd, we get 



kiT , kiT^ . 

6 = kitfj -h k 2 T -) y (0 + v) + -r (r ■+■ v) . 

00 (£ ^ d 



(23) 



3. Third order approximation in y 
In this case, 

y{t-T) = y-Ty + ]-T‘^y - \r^y^^ , 



and 



(24) 



8 = kill; + k 2 T + —y — -^(0 + u) -I- ~^—{r -1- u) — -^—{r 4- v) . (25) 

Xd Xd 2,Xd bxd 



4. First order approximation in 8 

If a time lag, T, exists on 8 instead of simply y, then 

8 = <^sat tanh 




where. 



8\ = <5o(t ~ T) = 8q — T8q . 
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This equation models a time lag associated with the entire application of the 
necessary corrective control action and not just its positional error. Using the 
above equations along with the linear equations 6 = 5\ and V’c = —y/^d, we 
get 

k k T' 

5 = kiij} + k^r H — k\Tr ^ — (•0 + v) — k 2 Tr . (26) 

^d 

5. First order approximation in both y and 5 

In a more general setting, we can assume that a time lag, Ti, exists in the 
control law S, and a different time lag, T 2 , is present in the processing of the 
positional error y. Assuming a first order approximation for both, we have 

~ Ti) — 60 — T 160 , 



and 



y{t - T 2 ) = y - T 2 y . 



Therefore, in this case using the above equations along with the linear equa- 
tions 6 and -^c = —y{i — T 2 )lxd ,we obtain 



k k T 

6 — ki'tl; -1- k 2 r -y — kiTir + v) — 

^d ^d 

k\T2 f V kiT\T2 f .V , ^ . 

-( 7 /; + I?) H (r + v) - k 2 TiT . 



Xd 



Xd 



(27) 



Using one of the above expressions, the linearized equations of motion can 
be written as 

Ip = r , 

V = dii^p Oi2V + dizr -h a\ 4 y -|- 616 , /2g\ 

r = 021 '0 "I" 0,22^ + 023^ + d24y "h > 

y = ‘lp + V , 

where all coefficients have been previously defined. 
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C. SERIES EXPANSIONS OF TIME LAG 

1. First order approximation in y 
In this case we have, 

y{t -T) = y - Ty . 



Substitution of (21) into (28) yields the following matrix system, 

0 



where. 



Tp 




■ 0 


0 


1 


V 




^21,1 


^22,1 


>l23,l 


r 




>131,1 


>132,1 


>133,1 


. y . 




1 


1 


0 



^21,1 — ttii + biki — 



0 



bikiT 





■ iP ■ 




V 




r 




. y . 



^22,1 = 0,12 — 



bikiT 

Xd 



^23,1 = ®13 + i>ik 2 

biki 



^24,1 — Ou + 



Xd 



^31,1 — 021 + ^2^1 ~ 



b 2 kiT 

Xd 



^32,1 = 022 — 



b 2 kiT 



^33,1 — ^23 + 62^2 

. , b 2 ki 

A34,1 = 024 H 

Xd 



2. Second order approximation in y 
In this case we have. 



1 . 



y{t-T) = y-Ty + ^T^y 



(29) 



17 



Substitution of (23) into (28) yields the following matrix system, 



where, 



■ i ) ■ 




■ 0 0 1 O' 




■ V’ ■ 


V 




■^21,2 -^22,2 ^23,2 ^24,2 




V 


f 




^31,2 ^32,2 ^33_2 ^34^2 




r 


. y . 




1 1 0 0 _ 




. y . 



^21,2 

^22,2 

^23,2 

■^24,2 

■^ 31,2 

^32,2 

^33,2 

■^ 34,2 



0-ii^d + bikjXd — bikjT 
Xd - O.SbikjT^ 
dj2^d ^ bjki'J^ 



xd - 0.56i/ciT2 
aizxd + bik2Xd + 0.55ifciT2 
Xd - O.bbikiT^ 
dl4Xd + ^1^1 



Xd - 0.5bikiT^ 



, , , b2k,T , b2k,T^ , 

<^21 + ® 2^1 1 n - di 21,2 



= 022 — 



Xd 2xrf 

b2kiT b2kiT^ 



+ 






2xd 



-■A22: 



, , b2kyT^ b2kiT^ ^ 

023 + ^2^2 d 1 — ^23,2 



— 024 + 



2xd 

62^1 b2k\T^ 



2Xn 



Xd 



+ 



2xd 



>^24,2 



3. Third order approximation in y 
A third order approximation in y is based on, 

y(t-T)=y-Ty + ir^jj - ■ 



( 30 ) 



Similar algebraic steps as in the previous two approximations result in the 
following eigenvalue problem. 



0 

0 

0 

0 

r-H 




■tp 


0 10 0 0 




V 


0 0 533^3 0 535^3 




r 


0 0 0 1 0 




y 


_ 0 0 5s3,3 0 555,3 . 




. ^2 . 
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( 31 ) 



■ 0 0 1 0 O' 




■ ■ 


0 0 0 0 1 




Vi 


^31,3 ^32,3 i 4 s 3^3 >134,3 ^435,3 




r 


110 0 0 




y 


. ^51,3 ^52,3 ^53,3 ^54,3 >^55, 3 . 




. ^’2 . 



where, v\ = u, V 2 is the rate of change of n, and the entries of the generalized 
eigenvalue problem (31) are given bellow. Higher order approximations in T 
can not produce any usable stability results, since the B matrix in (31) becomes 
singular. 



^31,3 


= 


Oil "b ^> 1^1 — 


bikiT 

Xd 


>^32, 3 


= 


6 ifcir 

Xd 




^33,3 


= 


ai3 + b\k2 + 


bikiT' 

2xd 


CO 

CO 




hiki 

Oi4 + 

Xd 


^35,3 


= 


hhT^ ^ 

2xd 




^51,3 


= 


021 + 62 /ci — 


b2kiT 

Xd 


> 152,3 


= 


62^16 

(^22 

Xd 




> 153,3 


= 


CL2Z + ^2^2 + 


b2kiT' 

2xd 


> 154,3 


= 


62 /ci 

024 + 

Xd 


^55,3 


= 


b2kiT‘^ 

2xd 




^33,3 


= 


b2kiT^ 

Qxd 




^35,3 


= 


633,3 




• 653,3 


= 


b2kiT^ 

6xd 
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•655,3 



b2kiT^ 



6xd 



4. First order approximation in 6 

Substitution of (26) into (28) yields the following matrix system, 

10 0 0 
0 1 623,4 0 
0 0 633,4 0 
0 0 0 1 





■ 1jj ■ 




■ 0 


0 


1 


0 




V 




-421,4 


-422,4 


-423,4 -424,4 




r 




-431,4 


-432,4 


-433,4 -434,4 




. y . 




1 


1 


0 


0 





■ ip ' 




V 




r 




. y . 



where, 



^21,4 — Oil + ~ 

bikiT 



bikiT 

Xd 



-422,4 = Ql2 - 






- 423,4 — ^^13 + b\k2 — bik\T 

b\ki 
Xd 

b2kiT 



-424,4 — Oi4 + 



-431,4 — <^21 + ^2^1 ~ 
^ b^kiT 

-^32,4 = ^22 — 



Xd 



Xd 



-433,4 = ^23 + ^2^2 ~ 62^16 

, , ^>2^1 

-^34,4 — ^^24 H 

^ d 

623.4 = b\k2T 

633.4 = 1 + ^2^26 

5. First order approximation in both y and 6 

Substitution of (27) into (28) yields the following matrix system. 



1 0 

0 622,5 
0 632,5 
0 0 



0 

623.5 

633.5 

0 



0 

0 

0 

1 





■ Ip ■ 






i] 






r 






. y . 





0 

-421,5 

-431,5 

1 



0 

-422,5 

-432,5 

1 



1 

-423,5 

-433,5 

0 



0 

-424,5 

-434,5 

0 



( 32 ) 





■ ip ■ 




V 




T 




. y . 



, (33) 
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where, 



^ 21,5 
^22,5 
^23,5 
>^24,5 
^31,5 
^ 32,5 
>^ 33, 5 
^ 34,5 
■ 522,5 

523.5 

532.5 
5 s 3,5 



h\kiT 2 

<^11 I ^1^1 

^ d 

bikiTi bikiT2 



(^12 ~ 






^d 



^13 + bik 2 ~ bikiTi + 



b\kiTiT 2 

^ d 



ai4 + 



feifc 



IM 

Xd 



®21 + ^2^1 ~ 



hkiTi b2kiT2 



0-22 — 



Xd Xd 

62^1X1 62^1X2 






a:<i 



®23 + & 2^2 ~ b2k\T\ + 

62^1 

^24 T 

Xd . 

^ _ felfciXiX2 

Xd 

bik 2 Ti 

62^1X1X2 

Xd 

1 + b 2 k 2 T\ 



62^1X1X2 

Xd 



D. EXACT COMPUTATION OF TIME LAG 

The previous analysis using Taylor series expansions of y{t — X) breaks 
down for approximations beyond third order, as the corresponding generalized 
eigenvalue problem becomes singular. Thus, in order to obtain an exact com- 
putation of the stability curves and check the validity of the calculations, a 
different technique will be performed in this section. This technique is based 
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on frequency response methods and it utilizes Nyquist’s criterion for stability 
[Friedland (1986)]. We write the system of equations (11) along with equations 
y = -0 + u and 6 = k\'il} + k 2 r + ~ 7"), in the Laplace domain, where 

y{t -T) — ^ ye“^* , (34) 

is the time delay operator. The characteristic equation of the system is, 

+ Q!i5^ + Q!2'5^ ”t" U 3 S + Q !4 + {02^^ 01^ d" 0o ) — ^ = 0 , (35) 

^ d 

where 

ai = — ai2 — 023 — hk2 , 

02 = — Oi 4 + 012023 — a22^\k2 — O 22 O 13 + Oi 2 & 2^2 ~ ^2^1 ~ O 21 , 

U3 = — O24 — 013U24 ~ 024^1^2 + O23O14 + 01462^2 ~ b\k\a22 ~ 

O11O22 d" ^2^1012 d- O21O12 , 

O4 = O21O14 + O12O24 — O22O14 — 011024 d- 62^1014 — 6ifci024 , 

/?2 = -6ifci , 

01 — ~ 01362^:1 + bik\a 2 z ~ t> 2 ki , 

00 — o.i2^2ki ~ bik\a22 ~ 52^1011 + 61^1021 , 

The characteristic equation is written as, 

1 H ^(^) ~ 0 ) 

Xd 

where 

(02S^ + 0\S 0(0e 

s'* + o;i 5 ^ + a2«^ d- 0:3s + 04 




(36) 
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is the open loop transfer function , and ^ denotes the effective position gain. 
For stability, we can utilize the Nyquist criterion which states that the critical 
value of Xd can be computed from, 



Xd' 



KGOw)! = 1 for axgG{ju) = -Tt 



(37) 



where j denotes the imaginary unit. The argument, argC?(ja;), of(36) is, 



a,rgG{jijj) = —u)T + tan 



-1 



f3iuj 



— tan 



a^oJ — aiu) 



Po — u)‘^ — + 04 

The set of equations (37) result in the critical visibility distance. 



(38) 



Xd = 









(39) 



{a^uj — aiw^)2 + (w^ — 02^^ + 04)^ 

The value of uj in(38) such that argG(ju>) = — tt is the value of the phase 
crossover frequency. After solving for the phase crossover frequency, the critical 
value of Xd is obtained by direct substitution of this value of w into equation 
(39). 



E. RESULTS AND DISCUSSION 

Typical results are presented in Figures 3 through 6. All results shown 
are nondimensional unless otherwise stated. The critical value of Xd versus ujn 
for zero time lag, and parametrized for different values of the damping ratio 
C is shown in Figure 3. Stability is ensured for values of Xd greater than its 
critical value. It can be seen that lower values of ujn require higher values 
of Xd for stability. This means that a less responsive helmsman will need a 
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1 





Figure 3: Critical xd versus for T 2 = 0 and various values of C, 




Figure 4: Critical xd versus u>n for C = 0.8 and various values of T 2 
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Figure 5: Critical xa versus iOn for Q — 0.8, T 2 — 5sec and various values of T\ 




Figure 6: Critical Xd versus for C = 0.8 and zero time lag: Channel effects 
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longer lookahead distance for a stable operation. Similar conclusions hold for 
variations in the damping ratio In this case, higher values of ^ correspond 
to better helmsman response which requires less lookahead distance. 

The effect of time lag T 2 is shown in Figure 4. Time lags are in seconds. It 
can be seen that reasonable amounts of time lag do not have a serious effect on 
stability, at least in a linear sense. Of course, an amount of time lag may have 
a serious effect on the transient response of the system as well as its ability to 
reject an external disturbance. This can only be established by a systematic 
series of numerical simulations. Similar conclusions hold for non-zero time 
lags Ti, as shown in Figure 5. 

Finally, the severe destabilizing effect of the canal is demonstrated by the 
results of Figure 6 where a comparison between canal and open water results is 
presented. It can be seen that an order of magnitude increase in the lookahead 
distance may be required if the same control parameters are to be used in both 
open water and a canal. 
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IV. NONLINEAR ANALYSIS 



A. INTRODUCTION 

It can be numerically confirmed that in all cases of stability loss of the pre- 
vious chapter, one pair of complex conjugate eigenvalues of the corresponding 
eigenvalue problem crosses transversally the imaginary axis. A situation like 
this in which a certain parameter is varied such that the real part of one pair 
of complex conjugate eigenvalues of the linearized system matrix crosses zero, 
results in the system leaving its steady state in an oscillatory manner. This 
loss of stability is called Hopf bifurcation and generically occurs in one of two 
ways, supercritical or subcritical. In the supercritical case, stable limit cycles 
are generated after the nominal straight line motion loses its stability. The 
amplitudes of these limit cycles are continuously increasing as the parameter 
distance from its critical value is increased. For small values of this criticality 
distance the resulting limit cycle is of small amplitude and differs little from 
the initial nominal state. In the subcritical ceise, however, stable limit cycles 
are generated before the nominal state loses its stability. Therefore, depend- 
ing on the initial conditions it is possible to diverge away from the nominal 
straight line path and converge towards a limit cycle even before the nominal 
motion loses its stability. This means that in the subcritical Hopf bifurcation 
case the domain of attraction of the nominal state is decreasing and in fact 
it shrinks to zero as the critical point is approached. Random external dis- 
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turbances of sufficient magnitude can throw the vehicle off to an oscillatory 
steady state even though the nominal state may still remain stable. After the 
nominal state becomes unstable, a discontinuous increase in the magnitude of 
motions is observed as there exist no simple stable nearby attractors for the 
vehicle to converge to. Distinction between these two qualitatively different 
types of bifurcation is, therefore, essential in design. The computational pro- 
cedure requires higher order approximations in the equations of motion and is 
the subject of the next section. 



B. DETAILED CALCULATIONS 



The nonlinear equations of motion are written as. 



r , 


(40) 


ail'0 -h Oi2^ "1“ d" 0,14V d" ) 


(41) 


a2\'lp d- 022V + 0,2ZT + a24V d- b2S' , 


(42) 


sin ip + V cos tjj , 


(43) 



where the control law is. 



6 = kiif -I- k 2 r + ki tan 



-iVit-T) 



Xd 



and, including saturation, 



6 — ^sat tanh. ( ^ ] . 



^sat , 



(44) 



(45) 



We perform a Taylor series expansion of the equations, keeping the first non- 
vanishing nonlinear terms. Due to port/starboard symmetry iri the problem 
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this means expansion to third order terms, 



(46) 

(47) 



sin^ = , cos'0 = 1 — , 




5 = kill) + k 2 r + —y{t -T)- - T) 

X d oX d 



(48) 



where for consistency, the term y{t — T) is approximated by its first order 



expansion in T, 

y{t - T) = y - Ty = y - Tip + \Tip^ - Tv + \Tvip^ . (49) 



Substitution of equations (46) to (49) into equations (40) to (45), produces 



the system. 



X = Ax + g(x) , 



(50) 



where the state variables vector is, 



X = [ijj.v.r.yY , 



A is the linearized matrix at equilibrium, and g(x) contains all third order 
terms. 

If T is the matrix of eigenvectors of A evaluated at the critical point Xd, 
the linear change of coordinates, 

x=Tz, z = T"^x, (51) 

transforms system (50) into its normal coordinate form, 

z = T-^ATz + T“^g(Tz) . (52) 
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At the Hopf bifurcation point, matrix T ^AT takes the form, 



T"^AT = 



0 —Wo 0 0 

Wo 0 0 0 

0 0 p 0 

0 0 0 g 






where wo is the imaginary part of the critical pair of eigenvalues, and the 
remaining two eigenvalues p and q are negative (or complex conjugate with 
negative real parts) . Therefore in the new system of coordinates (51) the 
dynamics of (50) are governed by a reduced two-dimensional system zi and 
Z 2 since the coordinates z^ and z^ correspond to the eigenvalues p and q and 
are asymptotically stable. For values of Xd close to the bifurcation point Xd, 
matrix T“^AT is. 



T”^AT = 



a'e 

Wo -t- w 
0 
0 



— (wq -I- W^£) 0 

a'e 0 

0 p + p'e 

0 0 



0 

0 

0 

Q + q'e 






where, e denotes the criticality difference. 



£ ^(icritical 



(53) 



w' = derivative of the real part of the critical eigenvalues with respect to e , 

a' = derivative of the imaginary part of the critical eigenvalues 

with respect to e , 

p — derivative of p with respect to e , 

q' = derivative of q with respect to e . 

Due to continuity, the eigenvalues p -f p'e and q + q'e remain negative for small 
nonzero values of e. Therefore, the coordinates 23, 24 correspond to negative 
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eigenvalues and are asymptotically stable. Center manifold theory predicts 
that the relationship between the critical coordinates zj, Z2 and the stable 
coordinates 23, 24 is at least of quadratic order. In fact, due to the symmetry 
in our problem the relationship is cubic, 

2i = 0 (zl 2I) , 22 = 0(2|, 2I) . 

It follows that because of this order, 23, 24 do not influence the third order 
Taylor series expansions, and can be dropped from the equations. Therefore, 
we can write the reduced system that describes the essential dynamics of ( 52 ) 
in the form, 

21 = a'ezi - (u>o + oj'e)z2 + ^ 1 (^ 1 , ^ 2 ) , (54) 

22 = (cjQ + + a'ez 2 + F 2 (zi, Z2) , ( 55 ) 

where Fi, F2 contain the third order terms, 

F\{z-^,Z2) = rn2^ + ri22i22 + ri 32 i 2 | + ri 42 ^ , (56) 

-^2(21, 22) = T2\zI + r22z\z2 + T2ZZ\zI + r 24^2 • (^ 7 ) 

The coefficients are computable from the previous Taylor expansions, and 
they are given below. 

The nonlinear expansion coefficients that are utilized in the definition 
of the cubic stability coefficient K, are given by the following; 

^11 = ^ 12^21 + ?^ 13^31 + ^ 14^41 + ’^ 12<^11 + ^ 13^^21 j 

ri2 = ^12^22 + n\z^Z2 + ^^14^42 + Wl2<^12 + ™13<^22 i 
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^13 — ^12^23 + ^13^33 + ^14^43 + ^12<^13 + ^13<^23 > 

ri4 = 7112^24 + ^13^34 + ^14^44 + ^12^^14 + ^13^24 t 

T2i = 7122^21 + ^23^31 + «24^41 + n22dll + 7l23<^21 , 

r22 = 7122^22 + ^23^32 + ^24^42 + ^22^12 + ^23^^22 j 

^23 = ^ 22^23 + ^23^33 + ^24^43 + ^22<^13 + ^23^23 j 

T24 = 7122^24 + 7123^34 + 7124^44 + 7l22<^14 + ^^23^^24 , 

where we denote, 

T = [rriij] , T"^ = [ 77 ^] , i,j = l,...,4. 



For numerical purposes the critical eigenvector of T must be normalized so 
that its first nonvanishing coefficient is identically equal to 1. The coefficients 
£ij are given by, 

^21 2 2 2 2 
, — ^77721 "1" ^^11^^31 ”b 

h 



2 2 2 2 
^ 41 ^'ipyy'^ll'^j^’i “1“ ^vrr^^2l^^3l 

^vyy^^2l'^ 4\ ”1“ ^rrj/ ^ryj/ 

~\~^^vr'^ll'^2l'^31 ^'ipvy'^ll'^21'^41 ”1“ ^^pry'^ll'^Al'^Sl 

3 3 3 3 

“l“(5^ry^^2l^^41^^3l ”1” “1” ^vvv'^21 S'rrr'^21 ^yyy'^41 ’ 

^ = S^^y(miiTn22 + 2777n7T7i27772l) + vv (777 i2m|i + 2777 11 777 1277722) 

0\ 

+ i5v,V>r(^11^32 + 2777 ii777 i27773i) + 6^rr{m\2m\-^ + 2777ii7773i77732) 
"l“^t/?t/>?/(777ii77742 "1“ 2777ii777 i27774l) (777^i777 i2 "1“ 2777ii7774i777 42) 

+^vvr(^21^32 "I" 2777317772177722) + ^tirr (7772277731 + 27773i777327772l) 
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+5yyy{m\-^m42 + 2miim2\m22) + + 2m4im42m2i) 

+^rry(^3lJ^42 + 2m4i77l3im32) + ^ryy (’Tl32?^4l 2m 4177142^131) 
+^V'i^r["7ll’^2l"732 + 77l3l(miim22 + 7771277121)] 

+^i/)uy[77llim2im42 + m4i(mnm22 + 77li2m2l)] 

+^^ry[777lim4im32 + m3i(mnm42 + 777i2m4i)] 

+^t,ry [777217774177732 + 77731 (m2im42 + 777227774i)j 

+^^^V,(3miimi2) + ^wu(3m2im22) + Srrr{‘ifn\imz2) + ^j/yy(3m4im42) 

^23 

— = 6^^„(mi27772i + 2miimi2m22) + ^^vv(777iim22 + 2mi2m2im22) 

0\ 

+<5^^r(777 i27773i + 2mumi2m32) + (5v-rr(777iim32 + 2mi2m3im32) 
+^^^^(777^277741 + 2miimi2m42) + ^V’yy("^42"^n + 2mi2m4im42) 

+ <5vvr(77722^31 2m2im2277732) + ^vTr{p^2\T>^‘z2 + 2m3im2277l32) 

+^„„j,(m 22777 41 + 2m2im2277742) + ^vyy (7772177742 + 2m4im42m22) 

+^rry(77732m4i + 2m42m3im32) + 5ryy(7773im42 + 2m4ir7742m32) 

+^V'ur[777 127772277731 + m32(mnm22 + 777i2m2l)] 

+^^uy[777 l2m227774i + m42(mnm22 + 777i2m2l)] 

+^V'r-y[777l277742m3i + m32(mnm42 + mi2m4i)] 

+5ury [777227774277731 + m32(m2im42 + m227774i)j 

+^^^^(37771177742) + ^wu( 377721777 22) + ^rrT’(37773l77732) + 

~~ — ^T^T^tf777 i277722 + 777 1277722 + ^T/j^r'777i277732 + *5T/jrr777 1277732 

0\ 

2 2 2 2 
22^^42 12^^42 22^^32 

2 2 2 2 
+^uuy7772277742 + *^uyy7772277742 + ^Try7773277742 + ^ryy 777 32 777 42 
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+<5^ur^ 12^22^32 ”1" ^i/’^’y^7ll2^22^42 + 12^42^X32 + 5t;^y^^22W^42^7^32 



^31 _ 

b 2 

^32 _ 

b 2 

^33 _ 

62 

^34 _ 

62 

£41 = 

^42 = 

^43 = 

£44 = 



•^‘2 "1“ byyyTJX^^ ^yyy^^42 * 

^21 
bi ’ 

^22 
61 ’ 

^23 

i>i ’ 

^24 

61 ’ 

--m^i (^m 2 i + -mii^ , 

/ 1 - 1 \ 
—mil ( 71^12^21 + ~7niim22 + ~7T^i2nT'ii ) j 

-mi2 (rniim22 + ^mi2m2i + ^mi2mn^ , 

-imi 2 (m 22 + 



The coefficients da are given by, 

1 



dll 

di 2 

di 3 

di 4 

d -21 

d 22 

d- 2 Z 

d 2 i 



Dv. 

J_ 

Dv, 

1 

Dvr 

1 

Dvr 

1 

Dvr 

1 

Dvr 

1 

Dvr 

1 

Dvr 



\c\\{Iz - Nf) + C 2 i(Tr - ma^c)] , 
•[ci2(-^z - A^r) + C22(^f ~ r^a^c)] , 
[ci3(-^z - A^r) + C23(^r ~ ma^c)] , 
[Cl4(-^z - AT^) + C24(^r “ raxc)] , 
[cii(7V^ - mxc) + C 2 i{m - y„)] , 
[cuiNv - mxc) + C 22 (m - Vi,)] , 
[ci3(A^i> - ruxo) + C 2 z{m - Vi)] , 
[ci4(7Vi, - mxa) + C24(m - Vi)] , 
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where 



Cii — 



+ 

Cl 2 = 
+ 
+ 

Cl3 = 
+ 
+ 

Ci4 = 
+ 

C 21 = 

+ 

C 22 = 
+ 
+ 

C 23 = 
+ 
+ 

C 24 = 

+ 



^ 3^2 ^2 

'X^vrr’^3'1^21 'Z^rvv^21^31 

0 2 2 

1 3 1 3 

+ g^yy^41 > 

iy„„^3m2im22 + ^V;r^(^7^3l”^22 + 2m3im32m2i) 

o z 

~I” 277131772-21^22) ”1” 'T^^‘0'03^11^12 

2 o 

g ■^y2/^^41^42 5 

^y„^^3m22r7721 + Wvrr{ml2m2l + 2771317773277722) 

o z 

(7772277731 + 2777327772177722) + “>^^,/,t/,3777 i2r77n 

^^2/2/2^42^41 ) 

O 

^3^2 ^2 

'^^ vvv '^22 

“^V’^t/>777i 2 + g^yyy77742 , 

~-^t;ui;7772i "f“ “■^^’"^^31^21 'Z^rvv'^21^^^ 

0 z z 

1 3 1 3 

— + gAryyy77741 , 

^A/'^^^37772i 77722 + ^A^i;rr (7773i77722 + 2777317773277721 ) 

0 Z 

-A^rt;t;(7772i77732 + 2777317772177722) + “A^i/;t/,t/»3777ii777i2 

1 2 

^-^2/2/2/ ^^41 ^42 5 
O 

-Nyyy3ml27n2l + -iV^rr(7773277721 + 2?77 3 177732777 22 ) 
^iV^t;t;(777227773i + 2777327772177722) + 7 3777 12777 n 

z o 

g A^2/2/2/^^42^41 3 

~^N yyyTfl22 “I" ~Z ^ vrr'^S2^'^'^ rvv'^22'^^*^ 

0 z z 

1 3 1 3 

g A/^t/ji/jt/> 77742 + g A^2/2/2/^42 * 
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The coefficients in a third order expansion of the control law are defined 



by 



^%l)yy = 

^VVT “ 
8 yTT ~ 
8yyy — 
8yyy = 

8 rry ~ 

6 ryy = 

8 %l)yr — 
8%jjyy — 

8 <ipry ~ 
8 yry — 



1 T 

72~(^i)^^ 2 + 0.5A:i 1 

St± Xd 



k^T^ 



■^sat 

1 



t -3 



T^ki 



-^k[{k',y+ , 

^sat 



1 



{k[yk 2 , 



<5sit 

1 



^2 

2 J 

^sat 



1 .kl , r/ci 



6 l 



sat 



1 



‘'1^2 
X J 






- 72 -(>= 2 ) k 2 , 



■'sat 

1 



^2 

2 f^ 2^2 ? 

^sat 
^sat 



^sat ^ 

1 ., 1 ^ Tk^ 

(2 ^2 2 + 

^sat 



^3 > 



Lfc2^ 

1 k"^ 

UJ^, 



5 l 



sat '^d 

1 



K 

f 1 f 



c2 ^^1^2^*^ 
8 sat 



1 .fci „r2fci 



^^1^2 2 3 , 
Xd x^ 



4. 

<^sat ^d 

“T^2^2^2— , 

^sat ^d 
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and 






6v 

6r 



+ 3x| 

1 1 



~ e s*^ 



^2 7 



sat 

1 



C _ ~ ~ "'1 _ 

'^yyy ~ -o _ •> 



1 A:f 



<5|at 3 



2fci 

3^’ 



3x^ 



k[ = 



ko = 



T 

ki - ki — , 
Xd 

i- ^ 

-ki — . 



C. NONLINEAR COEFFICIENT fC 

If we introduce polar coordinates in the form, 

z\ = R cos 6 , Z 2 = R sin 9 , (58) 

we can use (54) and (55) to produce an equation describing the rate of change 
of the radial coordinate R, 

R = a'eR + V{9)R^ , (59) 

and a similar equation in the rate of change of the angular coordinate 

9 = uo + uj'e + Q{e)R‘^ . (60) 

The system of equations (59) and (60) contains two variables, one of which, 
R, is slowly varying in time, whereas the other one, 0 is a fast variable. Then, 
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equation (59) can be averaged over one cycle in 9 to produce an equation with 
constant coefficients and similar stability properties, 

R = a'sR + ICR^, (61) 

where, 

1 

JC = — V{9)d9 = \ (3rn + ria + T22 + ^T2a) ■ (62) 

Z7T Jo 

Nontrivial equilibrium solutions of (61) correspond to limit cycles in the 
original system. The nontrivial equilibrium of (61), Rq, is given by, 

«o = ■ (63) 

In our problem, since the trivial equilibrium gains its stability for Xd > 
the coefficient a' is always negative. Therefore, it can be seen from (63) that 
a limit cycle will exist provided K and e have the same sign. The stability 
properties of this limit cycle can be determined by linearization of (61) around 
Rq. The linearized system is, 

R = -2a's{R - Ro) , (64) 

and its eigenvalue is, 

(3 = -2a’e . (65) 

We can see that the Floquet exponent (65) is negative if e is negative, which 
means that /C must be negative. Therefore, location and stability of limit 
cycles depends entirely on the cubic coefficient K. which is computable from 
(62). We can summarize our findings as. 
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omega n 



Figure 7: Nonlinear coefficient /C versus un for 5sat = 0.4 and various values of C 

• If /C < 0 then limit cycles exist for e < 0 (xa < and they are 

stable. 

• If /C > 0 then limit cycles exist for e > 0 (xd > and they are 

unstable. 



D. RESULTS AND DISCUSSION 

Typical results in terms of the nonlinear stability coefficient /C are shown in 
Figures 7 through 10. The general conclusion from this series of graphs is that 
the bifurcations are aU strongly subcritical. This means that any hnearized 
stability results should be viewed with extreme caution. Limit cycles exist 
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Figure 8: Nonlinear coefficient /C versus Uji for C = 0-8 and various values of <5gat 




omega n 

Figure 9: Nonlinear coefficient /C versus Un for C = 0-8, 6sat = 0.4, T\ = 0, and 
various values of (sec) 
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Figure 10: Nonlinear coefficient K. versus ojn with and without channel effects 

before stability in the linear sense is lost and a self sustained oscillation in the 
system may develop as a result of an external disturbance even if the nominal 
equilibrium state is still stable. The bifurcations become more subcritical; 
i.e., K. is more positive, for smaller values of C, as Figure 7 shows. Figure 8 
demonstrates the effect that the saturation limit 5gat has on the value of /C. 
Higher values of (Jgat, although are not related to the critical value of Xd result 
in significant changes in the value of JC. The general trend is that higher 5sat 
results in less subcritical bifurcations as evidenced by the overall decrease of 
JC. Figure 9 shows the effect of non-zero time lag on the nonlinear stability 
coefficient. It can be seen that, like the linear stabihty results, the effect is 
minimal. Finally, Figure 10 shows the canal effect on JC. This figure was 
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produced for zero time lag, C = 0.8, and Ssat = 0.4. It can be seen that the 
existence of a canal causes the bifurcations to be much more subcritical than 
the open water case. This demonstrates the severe destabilizing effect that the 
canal introduces in both the linearized and the nonlinear analysis. 
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V. CONCLUSIONS AND RECOMMENDATIONS 



This thesis presented a comprehensive study of linear and nonlinear stability 
properties of straightline motion of surface ships in confined waters. The 
classical maneuvering equations of motion incorporating canal suction effects, 
were coupled with appropriate navigation, gauidance, and control laws in order 
to mimic the helmsman’s behavior. The main conclusions from this study can 
be summarized as follows: 

1. There exists a critical preview distance for straightline positional stabil- 
ity. For values less than the critical distance, the system is unstable. 

2. Including canal effects, the critical preview distance may be an order of 
magnitude higher than in open water. If the same preview distance is to 
be used in both cases, the corresponding control law for canal maneuver- 
ing must be considerably more responsive than in open waters. 

3. The critical preview distance is monotonically decreasing for increasing 
control law responsiveness. This means that in order to accomodate 
smaller values of the preview distance, more responsive control laws are 
required . 

4. Physically realizable time lags do not seem to have a significant effect on 
the values of the preview distance. 

5. As the preview distance becomes less than its critical value, one pair of 
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complex conjugate eigenvalues of the linearized system matrix crosses 
the imaginary axis. This corresponds to a bifurcation to periodic so- 
lutions (Hopf bifurcation) and the system exhibits oscillatory behavior, 
also known as limit cycles. 

6. Higher order approximations in the equations of motion were utilized in 
order to assess the stability of the resulting limit cycles. It was found 
that in all cases, the limit cycles were unstable. This has the following 
implications: 

(a) It is possible for the system to lose its stability even before the critical 
preview distance is crossed. This means that all linearized stability 
analysis results should be viewed with extreme caution. 

(b) As the critical preview distance is crossed, it is expected that the 
system will develop limit cycles of large amplitude. This clearly 
presents a dangerous situation which should be avoided in practice, 
by appropriate changes in the design parameters. 

Recommendations for further research include the following: 

1. Simulation studies in order to verify the limit cycle behavior. 

2. Study of different ship characteristics and canal geometry. 
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APPENDIX 



The following is a list and description of the computer programs used in this 
thesis. The programs are written in FORTRAN. Complete printouts of the 
programs follow. Standard eigenvalue/eigenvector numerical analysis subrou- 
tines are required for all programs. 

• THESISl.FOR 

Calculation of critical Xd- First order approximation for time lag T 2 . 

• THESIS2.FOR 

Calculation of critical Xd- Second order approximation for time lag T 2 . 

• THESIS3.FOR 

Calculation of critical Xd- Third order approximation for time lag T 2 . 

. THESIS4.FOR 

Calculation of critical Xd- First order approximation for time lag T\. 

• THESIS5.FOR 

Calculation of critical Xd- First order approximation for time lags T\ and 

T2. 

• HOPF.FOR 

Calculation of the nonlinear cubic stability coefficient 1C. Requires the 

output of one of the previous programs as its input. 
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C PROGRAM THESISl.FOR (Time Delay-lst Order Approx. T2) 

C BIFURCATION ANALYSIS 

C 

PROGRAM THESIS 1 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DOUBLE PRECISION K1,K2,L,NR,NV,NDELTA,NPSI,NY,IZ,MASS, 
& NRDOT,NVDOT 

DIMENSION A(4,4) ,FV1(4) ,IV1(4) ,ZZZ(4,4) ,WR(4) ,VI(4) 

C 

OPEN (ll,FILE='BIFl.RESO 
OPEN (12,FILE=^BIF2.RESO 
OPEN (13,FILE=^BIF3.RESO 
open (15,f ile=^eig.resl O 
C 

C Vehicle Parameters : 

IZ =0.0 
L =528 
RHO =1.94 
XG =0.0 
MASS =0.0088 
U =1.0 
C 

YRDOT = 0.00000 
YVDOT =-0.00912 
YR =+0.00456 
YV =-0.01434 
YPSI = 0.01400 
YY = 0.02000 
YDELTA= 0.00278 
NRDOT =-0.00115 
NVDOT = 0.00000 
NR =-0.00296 
NV =-0.00460 
NPSI = 0.01000 
NY =-0.00250 
NDELTA=-0. 00139 
WRITE (*,1001) 

READ (*,*) WNMIN,WNMAX,IWN 
WRITE (*,1002) 

READ (*,*) XDMIN,XDMAX,IXD 
WRITE (*,1003) 

READ (*,*) ZETA 
WRITE(*,1100) 

READ (*,*) TL 
TL=TL*U/L 
C 

DVR =(IZ-NRDOT)*(MASS-YVDOT)- 
& (MASS*XG- YRDOT) * (MASS*XG-NVDOT) 
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AA 11= ( (IZ-NRDOT) *YPSI- (MASS*XG-YRDOT) *NPSI) /DVR 
AA12=( (IZ-NRDOT) *YV-(MASS*XG-YRDOT) *NV)/DVR 
AA21=( (NVDOT-MASS+XG) *YPSI+(MASS-YVDOT) *NPSI) /DVR 
AA22= ( (MASS-YVDOT) *NV- (MASS*XG-NVDOT) *YV) /DVR 
AA 13= ( ( IZ-NRDOT) * ( YR-MASS ) + 

& (YRDOT-MASS*XG) ♦ (NR-MASS+XG) ) /DVR 

AA23= ( (NVDOT-MASS*XG) * (YR-MASS) + 

& (MASS-YVDOT) * (NR-MASS*XG) ) /DVR 

AA14= ( (IZ-NRDOT) *YY+ (YRDOT-MASS*XG) *NY) /DVR 
AA24=( (NVDOT-MASS*XG) *YY+ (MASS-YVDOT) *NY) /DVR 
BBl =( (IZ-NRDOT) *YDELTA-(MASS*XG-YRDOT) *NDELTA) /DVR 
BB2 = ( (MASS-YVDOT) *NDELTA- (MASS*XG-NVDOT) *YDELTA) /DVR 

AN0M=AA23 
BN0M=BB2 
CN0M=AA21 
EPS =l.D-5 
ILMAX=1500 

DO 1 1=1, IWN 

WRITE (*,2001) I, IWN 

WN=WNMIN+ (I-l) * (WNMAX-WNMIN) / (IWN-1) 

Kl=- (WN**2 . DO/BNOM) - (CNOM/BNOM) 

K2=- ( ANOM+2 . DO*ZETA*WN) /BNOM 
DO 2 J=1,IXD 

XD=XDMIN+ ( J- 1 ) * (XDMAX-XDMIN) / ( IXD- 1 ) 

CALL LINEAR (K 1 , K2 , XD , AAl 1 , AA12 , AA13 , AA14 , AA21 , AA22 , AA23 , 
& AA24,BB1,BB2,A,TL) 

CALL RG(4,4,A,WR,WI,0,ZZZ,IV1,FV1,IERR) 

CALL DSTABL(DEOS,WR,WI,FREQ) 
write(15,*)DE0S,XD,WN 

IF (J.GT.l) GO TO 10 

DEOSOO=DEOS 

XDOO =XD 

LL=0 

GO TO 2 

10 DEOSNN=DEOS 
XDNN =XD 
PR=DEOSNN»DEOSOO 
IF (PR.GT.O.DO) GO TO 3 
LL=LL+1 

IF (LL.GT.3) STOP 1000 
IL=0 

XDO=XDOO 

XDN=XDNN 

DEOSO=DEOSOO 
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DEOSN=DEOSNN 
6 XDL=XDO 
XDR==XDN 
DEOSL=DEOSO 
DEOSR=DEOSN 
XD=(XDL+XDR)/2.D0 

CALL LINEAR (K1 , K2 , XD , AAl 1 , AA12 , AA13 , AA14 , AA21 , 
& AA22 , AA23 , AA24 , BBl , BB2 , A , TL) 

CALL RG(4,4,A,WR,WI,0,ZZZ,IV1,FV1,IERR) 

CALL DSTABL (DECS, WR,WI, FREQ) 

DEOSM=DEOS 

XDM=XD 

PRL=DEOSL*DEOSM 

PRR=DEOSR*DEOSM 

IF (PRL.GT.O.DO) GO TO 5 

XDO=XDL 

XDN=XDM 

DEOSO=DEOSL 

DEOSN=DEOSM 

IL=IL+1 

IF (IL.GT.ILMAX) STOP 3100 
DIF=DABS(XDL-XDM) 

IF (DIF.GT.EPS) GO TO 6 

XD=XDM 

GO TO 4 

5 IF (PRR.GT.O.DO) STOP 3200 
XDO=XDM 
XDN=XDR 
DEOSO=DEOSM 
DEOSN=DEOSR 
IL=IL+1 

IF (IL.GT.ILMAX) STOP 3100 
DIF=DABS(XDM-XDR) 

IF (DIF.GT.EPS) GO TO 6 
XD=XDM 

4 LLL=10+LL 

WRITE (LLL,*) XD,WN 
3 XDOO=XDNN 

DEOSOO=DEOSNN 
2 CONTINUE 
1 CONTINUE 



1001 


FORMAT 


(’ 


ENTER 


MIN, 


MAX, AND 


INCREMENTS 


OF 


WN 


(dimensionless) 


’) 


1002 


FORMAT 


(> 


ENTER 


MIN, 


MAX, AND 


INCREMENTS 


OF 


XD 


(dimensionless) 


’) 


1003 


FORMAT 


(' 


ENTER 


DAMPING RATION 


0 










1100 


FORMAT 


(’ 


ENTER 


TIME 


LAG TL (dimensional) 


0 








2001 


FORMAT 


(215) 
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END 



C 

SUBROUTINE DSTABL(DEOS,WR,WI .OMEGA) 

C 

C EVALUATES THE EIGENVALUE WITH THE MAXIMUM REAL PART 
C 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

DIMENSION WR(4),WI(4) 

DE0S=-1.0D+20 
DO 1 1=1,4 

IF (WR(I) .LT.DEOS) GO TO 1 
DEOS=WR(I) 

IJ=I 

1 CONTINUE 
OMEGA=WI(IJ) 

OMEGA=DABS (OMEGA) 

RETURN 

END 

C 

SUBROUTINE LINEAR(K1,K2,XD,AA11,AA12,AA13,AA14,AA21, 

& AA22.AA23,AA24,BB1,BB2,A.TL) 

C 

C FORMS THE LINEARIZED MATRIX A (time delay 1st order approximation 
C 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

DOUBLE PRECISION K1.K2 
DIMENSION A(4,4) 

A(l.l)=O.ODO 

A(1,2)=0.0D0 

A(1,3)=1.0D0 

A(1,4)=0.0D0 

A(2 , 1) =AA11+BB1*K1-BB1*K1*TL/XD 

A(2,2)=AA12-BB1*K1*TL/XD 

A(2,3)=AA13+BB1*K2 

A(2 , 4) =AA14+BB1*K1/XD 

A(3,1)=AA21+BB2*K1-BB2*K1*TL/XD 

A(3.2)=AA22-BB2*K1*TL/XD 

A(3,3)=AA23+BB2*K2 

A (3 , 4) =AA24+BB2*K1/XD 

A(4,1)=1.0D0 

A(4,2)=1.0D0 

A(4,3)=0.0D0 

A(4.4)=0.0D0 

RETURN 

END 
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c 

c 

c 



PROGRAM THESIS2.F0R (Time Delay-2nd Order Approx T2) 
BIFURCATION ANALYSIS 

PROGRAM THESIS2 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DOUBLE PRECISION K1 ,K2,L ,NR,NV,NDELTA ,NPSI , NY, IZ ,MASS , 
& NRDOT,NVDOT 

DIMENSION A(4,4) ,FV1 (4) , IVl (4) ,ZZZ(4,4) ,WR(4) ,WI(4) 



C 

C 



OPEN 


dl,FILE=’BIFl.RES’) 


OPEN 


(12,FILE='BIF2.RES') 


OPEN 


(13,FILE=’BIF3.RESO 


open 


(15,file=’eig-resl’) 


Vehicle Paxameters: 


IZ 


=0.0 


L 


=528 


RHO 


=1.94 


XG 


=0.0 


MASS 


=0 . 0088 


U 


=1.0 


YRDOT 


= 0.00000 


YVDOT 


=-0.00912 


YR 


=+0.00456 


YV 


=-0.01434 


YPSI 


= 0.01400 


YY 


= 0.02000 


YDELTA= 0.00278 


NRDOT 


=-0.00115 


NVDOT 


= 0.00000 


NR 


=-0.00296 


NV 


=-0.00460 


NPSI 


= 0.01000 


NY 


=-0.00250 


NDELTA=“0. 00139 


WRITE 


(♦,1001) 


READ 


(♦,♦) WNMIN.WNMAX 


WRITE 


(♦,1002) 


READ 


(♦,♦) XDMIN,XDMAX 



WRITE (*,1003) 
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ZETA 



READ (*,*) 

WRITE(*,1100) 

READ (*,*) TL 
TL=TL*U/L 
C 

DVR =(IZ-NRDOT)*(MASS-YVDOT)- 
& (MASS*XG-YRDOT)*(MASS*XG-NVDOT) 

AAl 1= ( ( IZ-NRDOT) *YPSI- (MASS*XG-YRDOT) *NPS I ) /DVR 
AA12=( (IZ-NRDOT)*YV-(MASS*XG-YRDOT) *NV)/DVR 
AA21= ( (NVDOT-MASS*XG) *YPSI+(MASS-YVDOT) *NPSI)/DVR 
AA22= ( (MASS-YVDOT) *NV- (MASS*XG-NVDOT) *YV) /DVR 
AA13=( (IZ-NRDOT) ♦(YR-MASS)+ 
k (YRDOT-MASS*XG) ♦ (NR-MASS*XG) ) /DVR 

AA23=( (NVDOT-MASS*XG) * (YR-MASS) + 
k (MASS-YVDOT) *(NR-MASS*XG)) /DVR 

AA14=( (IZ-NRDOT) *YY+(YRDOT-MASS*XG) *NY) /DVR 
AA24= ( (NVDOT-MASS*XG) *YY+ (MASS-YVDOT) *NY) /DVR 

BBl =( (IZ-NRDOT) *YDELTA- (MASS*XG-YRDOT) *NDELTA) /DVR 
BB2 = ( (MASS-YVDOT) *NDELTA- (MASS*XG-NVDOT) *YDELTA) /DVR 
C 

AN0M=AA23 

BN0M=BB2 

CN0M=AA21 

EPS =l.D-5 
ILMAX=1500 
C 

DO 1 1=1, IWN 

WRITE (*,2001) I, IWN 
WN=WNMIN+(I-1) * (WNMAX-WNMIN)/ (IWN-1) 

Kl=- (WN**2 .DO/BNOM) -(CNOM/BNOM) 

K2=- (ANOM+2 . DO*ZETA*WN) /BNOM 

DO 2 J=1,IXD 

XD=XDMIN+ ( J- 1 ) * (XDMAX-XDMIN) / (IXD- 1 ) 

C 

CALL LINEAR (K 1 , K2 , XD , AAl 1 , AA12 , AA13 , AA14 , AA2 1 , AA22 , AA23 , 
k AA24,BB1,BB2,A,TL) 

CALL RG(4,4,A,WR,WI,0,ZZZ,IV1,FV1,IERR) 

CALL DSTABL(DEOS,WR,WI,FREQ) 
write ( 15 , ♦ ) DEOS , XD , WN 
C 

IF (J.GT.l) GO TO 10 

DEOSOO=DEOS 

XDOO =XD 
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LL=0 
GO TO 2 

10 DEOSNN=DEOS 
XDNN =XD 
PR=DEOSNN*DEOSOO 
IF (PR.GT.O.DO) GO TO 3 
LL=LL+1 

IF (LL.GT.3) STOP 1000 
IL=0 

XDO=XDOO 
XDN=XDNN 
DEOSO=DEOSOO 
DEOSN=DEOSNN 
6 XDL=XDO 

XDR=XDN 
DEOSL=DEOSO 
DEOSR=DEOSN 
XD=(XDL+XDR)/2.D0 
C 

CALL LINEAR (K1 , K2 , XD , AAl 1 , AA12 , AA13 , AA14 , AA2 1 , 
& AA22 , AA23 , AA24 , BBl , BB2 , A , TL) 

CALL RG(4,4,A,WR,WI,0,ZZZ,IV1,FV1,IERR) 

CALL DSTABL(DEOS,WR,WI,FREQ) 

DEOSM=DEOS 

XDM=XD 

PRL=DEOSL*DEOSM 

PRR=DEOSR*DEOSM 

IF (PRL.GT.O.DO) GO TO 5 

XDO=XDL 

XDN=XDM 

DEOSO=DEOSL 

DEOSN=DEOSM 

IL=IL+1 

IF (IL.GT.ILMAX) STOP 3100 
DIF=DABS(XDL-XDM) 

IF (DIF.GT.EPS) GO TO 6 

XD=XDM 

GO TO 4 

5 IF (PRR.GT.O.DO) STOP 3200 
XDO=XDM 
XDN=XDR 
DEOSO=DEOSM 
DEOSN=DEOSR 
IL=IL+1 

IF (IL.GT.ILMAX) STOP 3100 
DIF=DABS(XDM-XDR) 

IF (DIF.GT.EPS) GO TO 6 
XD=XDM 
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4 



LLL=10+LL 

WRITE (LLL,*) XD,WN 
3 XDOO=XDNN 

DEOSOO=DEOSNN 
2 CONTINUE 
1 CONTINUE 
C 



1001 


FORMAT 


(’ 


ENTER 


MIN, 


MAX, AND INCREMENTS 


OF 


WN 


(dimensionless) 0 


1002 


FORMAT 


(’ 


ENTER 


MIN, 


MAX, AND INCREMENTS 


OF 


XD 


(dimensionless) ’) 


1003 


FORMAT 


(’ 


ENTER 


DAMPING RATIO’) 








1100 


FORMAT 


(’ 


ENTER 


TIME 


LAG TL (dimensional) 


0 






2001 


FORMAT 


(215) 













END 

C 

SUBROUTINE DSTABL(DEOS,WR,WI .OMEGA) 

C 

C EVALUATES THE EIGENVALUE WITH THE MAXIMUM REAL PART 
C 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

DIMENSION WR(4) ,WI(4) 

DE0S=-1.0D+20 
DO 1 1=1,4 

IF (WR(I) .LT.DEOS) GO TO 1 
DEOS=WR(I) 

IJ=I 

1 CONTINUE 
OMEGA=WI(IJ) 

OMEGA=DABS (OMEGA) 

RETURN 

END 

C 

SUBROUTINE LINEAR(K1 ,K2,XD, AAll ,AA12,AA13, AA14.AA21 , 

& AA22,AA23,AA24,BB1,BB2,A,TL) 

C 

C FORMS THE LINEARIZED MATRIX A (time delay 1st order approximation) 
C 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

DOUBLE PRECISION K1.K2 
DIMENSION A (4, 4) 

A(1,1)=O.ODO 

A(1,2)=O.ODO 

A(1,3)=1.0D0 

A(1,4)=0.0D0 

A (2 , 1) = (AAl 1 *XD+BB1*K1*XD-BB1*K1*TL) / (XD-0 . 50D0*BB1*K1*TL*TL) 
A(2,2)=(AA12*XD-BB1*K1*TL)/(XD-0.50D0*BB1*K1*TL*TL) 

A (2 , 3) = (AAl 3*XD+BBl*K2*XD+0 . 50D0*BB1 *K1*TL*TL) / 

& (XD-0.50D0*BB1*K1*TL*TL) 

A(2,4)=(AA14*XD+BB1*K1)/(XD-0.50D0*BB1*K1*TL*TL) 
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A (3 , 1) =AA21+BB2*K1-BB2*K1*TL/XD+ (BB2*K1*TL»TL/ (2 . ODO»XD) ) *A (2 , 1) 
A (3 , 2) =AA22-BB2*K1*TL/XD+ (BB2*K1*TL*TL/ ( 2 . ODO*XD) ) *A (2 , 2) 

A (3 , 3) =AA23+BB2*K2+ (BB2»K1*TL*TL/ (2 . ODO*XD) )+ 

& (BB2*K1*TL*TL/(2.0*XD))* A (2, 3) 

A (3 , 4) =AA24+BB2*K1/XD+ (BB2*K1 *TL*TL/ (2 . ODO*XD) ) * A ( 2 . 4) 

A(4,1)=1.0D0 

A(4,2)=1.0D0 

A(4,3)=0.0D0 

A(4,4)=0.0D0 

RETURN 

END 
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C PROGRAM THESIS3.F0R (Time Delay-3rd Order Approx T2) 

C BIFURCATION ANALYSIS 

C 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

DOUBLE PRECISION K1 ,K2,L,NR,NV,IZ,MASS ,NDELTA,NPSI ,NY, 
& NRDOT.NVDOT 

DIMENSION A(5,5) ,B(5,5) ,FV1(5) ,IV1(5) ,ZZZ(5,5) ,ALFR(5) , 
& ALFI(5),BETA(5) ,WR(5) ,WI(5) 

C 

OPEN (11,FILE=’BIF1.RES’) 

OPEN (12,FILE=’BIF2.RES’) 

OPEN (13,FILE='’BIF3.RES’) 

C 

C Vehicle Parameters: 

IZ =0.0 
L =528 
RHO =1.94 
XG =0.0 
MASS =0.0088 
U =1.0 
C 

YRDOT = 0.00000 
YVDOT =-0.00912 
YR =+0.00456 
YV =-0.01434 
YPSI = 0.01400 
YY = 0.02000 
YDELTA= 0.00278 
NRDOT =-0.00115 
NVDOT = 0.00000 
NR =-0.00296 
NV =-0.00460 
NPSI = 0.01000 
NY =-0.00250 
NDELTA=-0. 00139 
WRITE (*,1001) 

READ (*,*) WNMIN.WNMAX.IWN 
WRITE (*,1002) 

READ (*,*) XDMIN,XDMAX,IXD 
WRITE (*,1003) 

READ (*,*) ZETA 
WRITE(*,1100) 

READ (*,*) TL 
TL=TL*U/L 
C 

DVR =(IZ-NRDOT)*(MASS-YVDOT)- 
& (MASS*XG-YRDOT)*(MASS*XG-NVDOT) 
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AAl 1= ( ( IZ-NRDOT) *YPSI- (MASS*XG-YRDOT) *NPSI ) /DVR 
AA12= ( ( IZ-NRDOT) *YV- (MASS+XG-YRDOT) ♦NV) /DVR 
AA21= ( (NVDOT-MASS*XG) *YPSI+ (MASS-YVDOT) *NPSI) /DVR 
AA22= ( (MASS-YVDOT) +NV- (MASS*XG-NVDOT) *YV) /DVR 
AA13= ( ( IZ-NRDOT) * ( YR-MASS) + 

& ( YRDOT-MASS*XG) ♦ (NR-MASS*XG) ) /DVR 

AA23=( (NVDOT-MASS+XG) * (YR-MASS) + 

& (MASS-YVDOT) *(NR-MASS*XG)) /DVR 

AA14=( (IZ-NRDOT) *YY+(YRDOT-MASS*XG)*NY)/DVR 
AA24= ( (NVDOT-MASS*XG) *YY+ (MASS-YVDOT) *NY) /DVR 

BBl = ( (IZ-NRDOT) *YDELTA- (MASS*XG-YRDOT) *NDELTA) /DVR 
BB2 = ( (MASS-YVDOT) *NDELTA- (MASS*XG-NVDOT) *YDELTA) /DVR 

AN0M=AA23 

BN0M=BB2 

CN0M=AA21 

EPS =l.D-5 
ILMAX=1500 

DO 1 1=1, IWN 

WRITE (*,2001) I, IWN 

WN=WNMIN+ ( I- 1) * ( WNMAX-WNMIN) / ( IWN- 1 ) 

Kl=- (WN**2 . DO/BNOM) - (CNOM/BNOM) 

K2=- ( ANOM+2 . DO *ZETA*WN) /BNOM 
DO 2 J=1,IXD 

XD=XDMIN+ ( J- 1 ) * (XDMAX-XDMIN ) / ( IXD- 1 ) 

CALL LINEAR(K1,K2,XD,AA11,AA12,AA13,AA14, 

AA2 1 , A A22 , A A23 , A A24 , BB 1 , BB2 , A , B , TL ) 

CALL RGG(5,5,A,B,ALFR,ALFI,BETA,0,ZZZ,IERR) 

DO 11 IJE=1,5 

WR ( IJE) =ALFR ( IJE) /BETA ( IJE) 

WI (IJE)=ALFI (IJE) /BETA (IJE) 

CONTINUE 

CALL DSTABL(DEOS,WR,WI,FREQ) 

IF (J.GT.l) GO TO 10 
DEOSOO=DEOS 
XDOO =XD 
LL=0 
GO TO 2 
DEOSNN=DEOS 
XDNN =XD 
PR=DEOSNN*DEOSOO 
IF (PR.GT.O.DO) GO TO 3 
LL=LL+1 



C 

& 

11 

C 
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IF (LL.GT.3) STOP 1000 
IL=0 

XDO=XDOO 
XDN=XDNN 
DEOSO=DEOSOO 
DEOSN=DEOSNN 
6 XDL=XDO 
XDR=XDN 
DEOSL=DEOSO 
DEOSR=DEOSN 
XD=(XDL+XDR)/2.D0 

CALL LINEAR (K1 , K2 , XD , AAl 1 , AA12 , AA13 , AA14 , AA2 1 , AA22 , 
k AA23,AA24,BB1,BB2,A,B,TL) 

CALL RGG (5 , 5 , A , B , ALFR , ALFI , BETA , 0 , ZZZ , lERR) 

DO 12 IJE=1,5 

WR ( IJE) =ALFR (IJE) /BETA ( IJE) 

WI (IJE)=ALFI (IJE) /BETA (IJE) 

12 CONTINUE 

CALL DSTABL(DEOS,WR,WI,FREQ) 

DEOSM=DEOS 

XDM=XD 

PRL=DEOSL*DEOSM 

PRR=DEOSR*DEOSM 

IF (PRL.GT.O.DO) GO TO 5 

XDO=XDL 

XDN=XDM 

DEOSO=DEOSL 

DEOSN=DEOSM 

IL=IL+1 

IF (IL.GT.ILMAX) STOP 3100 
DIF=DABS(XDL-XDM) 

IF (DIF.GT.EPS) GO TO 6 

XD=XDM 

GO TO 4 

5 IF (PRR.GT.O.DO) STOP 3200 

XDO=XDM 
XDN=XDR 
DEOSO=DEOSM 
DEOSN=DEOSR 
IL=IL+1 

IF (IL.GT.ILMAX) STOP 3100 
DIF=DABS(XDM-XDR) 

IF (DIF.GT.EPS) GO TO 6 
XD=XDM 

4 LLL=10+LL 

WRITE (LLL,*) XD,WN 
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3 XDOO=XDNN 

DEOSOO=DEOSNN 
2 CONTINUE 
1 CONTINUE 



C 



1001 

1002 

1003 

1100 

2001 



FORMAT (’ ENTER MIN, MAX, AND INCREMENTS OF WN (dimensionless)’) 
FORMAT (’ ENTER MIN, MAX, AND INCREMENTS OF XD (dimensionless)’) 
FORMAT (’ ENTER DAMPING RATIO’) 

FORMAT (’ ENTER TIME LAG TL (dimensional)’) 

FORMAT (215) 

END 



SUBROUTINE DSTABL(DEOS,WR,WI, OMEGA) 
IMPLICIT DOUBLE PRECISION (A-H,0-Z) 
DIMENSION WR(5),WI(5) 

DE0S=-1.0D+20 
DO 1 1=1,5 

IF (WR(I) .LT.DEOS) GO TO 1 
DEOS=WR(I) 

U=I 



1 CONTINUE 
OMEGA=WI(IJ) 
OMEGA=DABS (OMEGA) 
RETURN 
END 



SUBROUTINE LINEAR (K1 , K2 , XD , AAl 1 , AA12 , AA13 , AA14 , AA21 , AA22 , AA23 , AA24 
& ,BB1,BB2,A,B,TL) 



IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DOUBLE PRECISION K1,K2 
DIMENSION A(5,5) ,B(5,5) 

A(1,1)=0.0D0 

A(1,2)=0.0D0 

A(1,3)=1.0D0 

A(1,4)=0.0D0 

A(1,5)=0.0D0 

A(2,1)=0.0D0 

A(2,2)=0.0D0 

A(2,3)=0.0D0 

A(2,4)=0.0D0 

A(2,5)=1.0D0 

A(3,1)=AA11-(BB1*K1*TL/XD)+BB1*K1 

A(3,2)=AA12-BB1*K1*TL/XD 

A (3 , 3)=AA13+BB1*K2+ (BB1*K1*TL*TL/ (2 . DO*XD) ) 
A(3,4)=AA14+BB1*K1/XD 
A(3,5)=-1.D0+(BB1*K1*TL*TL/(2.D0*XD)) 
A(4,1)=1.0D0 
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A(4,2)=1.0D0 

A(4,3)=0.0D0 

A(4,4)=0.0D0 

A(4,5)=0.0D0 

A (5 , 1) =AA21+BB2*K1 -BB2*K1*TL/XD 
A(5 , 2) =AA22-BB2*K1*TL/XD 

A(5 , 3) =AA23+BB2*K2+ (BB2*K1*TL*TL/ (2 . DO*XD) ) 

A(5,4)=AA24+BB2*K1/XD 

A(5 , 5) = (BB2*K1*TL*TL/ (2 . DO+XD) ) 

C 

B(l,l)=l-ODO 

B(1,2)=0.0D0 

B(1,3)=0.0D0 

B(1,4)=0.0D0 

B(1,5)=0.0D0 

B(2,1)=0.0D0 

B(2,2)=1.0D0 

B(2,3)=0.0D0 

B(2,4)=0.0D0 

B(2,5)=0.0D0 

B(3,1)=0.0D0 

B(3,2)=0.0D0 

B (3 , 3) = (BBl *K1*TL*TL*TL/ (6 . DO*XD) ) 
B(3,4)=0.0D0 

B (3 , 5) = (BBl *K1*TL*TL*TL/ (6 . DO*XD) ) 

B(4,1)=0.0D0 

B(4,2)=0.0D0 

B(4,3)=0.0D0 

B(4,4)=1.0D0 

B(4,5)=0.0D0 

B(5,1)=0.0D0 

B(5,2)=0.0D0 

B (5 , 3) =1 . ODO+ (BB2*K1*TL*TL*TL/ (6 . DO*XD) ) 
B(5,4)=0.0D0 

B (5 , 5) = (BB2*K1*TL*TL*TL/ (6 . DO*XD) ) 

RETURN 

END 
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C PROGRAM THESIS4.FOR (Time Delay-lst Order Approx in delta ,T1) 
C BIFURCATION ANALYSIS 
C 



IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DOUBLE PRECISION K1,K2,L,NR,NV,IZ,MASS,NDELTA,NPSI.NY. 

& NRDOT.NVDOT 

DIMENSION A(4,4) ,B(4,4) ,FV1(4) ,IV1(4) ,ZZZ(4,4) ,ALFR(4) . 
& ALFK4) ,BETA(4) ,WR(4) ,WI(4) 



OPEN (11,FILE=’BIF1.RES’) 
OPEN (12,FILE=’BIF2.RES’) 
OPEN (13,FILE='BIF3.RES’) 
C 

C Vehicle Parameters: 



IZ 


=0.0 


L 


=528 


RHO 


=1.94 


G 


=32.2 


XG 


=0.0 


MASS 


=0.0088 


U 


=3.0 



YRDOT 


= 0.00000 


YVDOT 


=-0.00912 


YR 


=+0.00456 


YV 


=-0.01434 


YPSI 


= 0.01400 


YY 


= 0.02000 



YDELTA= 0.00278 
NRDOT =-0.00115 
NVDOT = 0.00000 
NR =-0.00296 
NV =-0.00460 

NPSI = 0.01000 

NY = -0.00250 

NDELTA=-0 . 00139 

WRITE (*,1001) 

READ (*,*) WNMIN.WNMAX.IWN 

WRITE (*,1002) 

READ (*,*) XDMIN,XDMAX,IXD 

WRITE (*,1003) 
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ZETA 



READ (*,*) 

WRITE(*,1100) 

READ (*,*) TL 
TL=TL*U/L 

DVR =(IZ-NRDOT)*(MASS-YVDOT)- 
& (MASS*XG-YRDOT)*(MASS*XG-NVDOT) 

AA11=((IZ-NRD0T)*YPSI-(MASS*XG-YRD0T)*NPSI)/DVR 
AA12= ( (IZ-NRDOT) *YV- (MASS*XG-YRDOT) *NV) /DVR 
AA21=((NVD0T-MASS*XG)*YPSI+(MASS-YVD0T)*NPSI)/DVR 
AA22=( (MASS-YVDOT) *NV- (MASS*XG-NVDOT) *YV) /DVR 
AA13= ( ( IZ-NRDOT) * ( YR-MASS) + 

& (YRDOT-MASS+XG) * (NR-MASS*XG) ) /DVR 

AA23=( (NVDOT-MASS*XG) * (YR-MASS) + 

& (MASS-YVDOT) *(NR-MASS*XG)) /DVR 

AA14=( (IZ-NRDOT) *YY+(YRDOT-MASS*XG)*NY)/DVR 
AA24=( (NVDOT-MASS+XG) *YY+ (MASS-YVDOT) *NY) /DVR 

BBl =( (IZ-NRDOT) *YDELTA- (MASS *XG-YRDOT)*NDELTA) /DVR 
BB2 =( (MASS-YVDOT) *NDELTA- (MASS*XG-NVDOT) *YDELTA) /DVR 

AN0M=AA23 

BN0M=BB2 

EPS =l.D-5 
ILMAX=1500 

DO 1 1=1, IWN 

WRITE (*.2001) I, IWN 

WN=WNMIN+ ( I- 1 ) * ( WNMAX-WNMIN ) / ( IWN- 1 ) 

K1=-WN**2.D0/BN0M 

K2=- ( ANOM+2 . DO*ZETA*WN) /BNOM 

DO 2 J=1,IXD 

XD=XDMIN+ ( J- 1) * (XDMAX-XDMIN) / ( IXD- 1) 

CALL LINEAR(K1 , K2 ,XD , AAl 1 , AA12 , AA13 , AA14 , 

AA2 1 . A A22 , AA23 , AA24 . BB 1 , BB2 , A , B . TL) 

CALL RGG (4 , 4 , A , B , ALFR , ALFI , BETA , 0 , ZZZ , lERR) 

DO 11 IJE=1,4 

WR ( IJE) =ALFR ( I JE) /BETA ( I JE) 

WI ( IJE) =ALFI (IJE) /BETA ( IJE) 

CONTINUE 

CALL DSTABL(DEOS,WR,WI,FREQ) 

IF (J.GT.l) GO TO 10 
DEOSOO=DEOS 
XDOO =XD 
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LL=0 
GO TO 2 

10 DEOSNN=DEOS 
XDNN =XD 
PR=DEOSNN*DEOSOO 
IF (PR.GT.O.DO) GO TO 3 
LL=LL+1 

IF (LL.GT.3) STOP 1000 
IL=0 

XDO=XDOO 
XDN=XDNN 
DEOSO=DEOSOO 
DEOSN=DEOSNN 
6 XDL=XDO 
XDR=XDN 
DEOSL=DEOSO 
DEOSR=DEOSN 
XD=(XDL+XDR)/2.D0 

CALL LINEAR(K1,K2,XD,AA11,AA12,AA13,AA14,AA21,AA22, 
& AA23,AA24,BB1,BB2,A,B,TL) 

CALL RGG (4 , 4 , A , B , ALFR , ALFI , BETA , 0 , ZZZ , lERR) 

DO 12 IJE=1,4 

WR ( IJE) =ALFR (I JE) /BETA ( IJE) 

WI (IJE) =ALFI (IJE) /BETA (IJE) 

12 CONTINUE 

CALL DSTABL(DEOS,WR,WI,FREQ) 

DEOSM=DEOS 

XDM=XD 

PRL=DEOSL*DEOSM 

PRR=DEOSR*DEOSM 

IF (PRL.GT.O.DO) GO TO 5 

XDO=XDL 

XDN=XDM 

DEOSO=DEOSL 

DEOSN=DEOSM 

IL=IL+1 

IF (IL.GT.ILMAX) STOP 3100 
DIF=DABS(XDL-XDM) 

IF (DIF.GT.EPS) GO TO 6 

XD=XDM 

GO TO 4 

5 IF (PRR.GT.O.DO) STOP 3200 
XDO=XDM 
XDN=XDR 
DEOSO=DEOSM 
DEOSN=DEOSR 
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IL=IL+1 

IF (IL.GT.ILMAX) STOP 3100 
DIF=DABS(XDM-XDR) 

IF (DIF.GT.EPS) GO TO 6 
XD=XDM 



4 LLL=10+LL 

WRITE (LLL,*) XD.WN 
3 XDOO=XDNN 

DEOSOO=DEOSNN 
2 CONTINUE 
1 CONTINUE 



C 



1001 

1002 

1003 

1100 

2001 



FORMAT (’ ENTER MIN, MAX, AND INCREMENTS OF WN (dimensionless)’) 
FORMAT (’ ENTER MIN, MAX, AND INCREMENTS OF XD (dimensionless)’) 
FORMAT (’ ENTER DAMPING RATIO’) 

FORMAT (’ ENTER TIME LAG TL (dimensional)’) 

FORMAT (215) 

END 



SUBROUTINE DSTABL(DEOS,WR,WI , OMEGA) 
IMPLICIT DOUBLE PRECISION (A-H,0-Z) 
DIMENSION WR(4),WI(4) 

DE0S=-1.0D+20 



DO 1 1=1,4 

IF (WR(I).LT.DEOS) GO TO 1 
DEOS=WR(I) 

U=I 



1 CONTINUE 
OMEGA=WI(IJ) 
OMEGA=DABS (OMEGA) 
RETURN 
END 



SUBROUTINE LINEAR (K 1 , K2 , XD , AA 1 1 , AA12 , AA13 , AA14 , AA2 1 , AA22 , AA23 , A A24 
k ,BB1,BB2,A,B,TL) 



IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DOUBLE PRECISION K1,K2 

DIMENSION A(4,4) ,B(4,4) 

A(1,1)=O.ODO 

A(1,2)=0.0D0 

A(1,3)=1.0D0 

A(1,4)=O.ODO 



A(2,1)=AA11+BB1*K1-BB1*K1*TL/XD 
A (2 , 2) =AA12-BB1*K1*TL/XD 
A(2,3)=AA13+BB1*K2-BB1*K1*TL 
A(2,4)=AA14+BB1*K1/XD 
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A(3, 1)=AA21+BB2+K1-BB2*K1*TL/XD 
A (3 , 2) =AA22-BB2*K1*TL/XD 
A(3,3)=AA23+BB2*K2-BB2+K1*TL 
A (3 , 4) =AA24+BB2*K1/XD 

A(4,1)=1.0D0 

A(4,2)=1.0D0 

A(4,3)=0.0D0 

A(4,4)=0.0D0 

B(1,1)=1.0D0 

B(1,2)=0.0D0 

B(1,3)=0.0D0 

B(1,4)=0.0D0 

B(2,1)=0.0D0 

B(2,2)=1.0D0 

B(2,3)=BB1*K2+TL 

B(2,4)=0.0D0 

B(3,1)=0.0D0 
B(3,2)=0.0D0 
B (3 , 3) =1 . 0D0+ (BB2*K2*TL) 
B(3,4)=0.0D0 

B(4,1)=0.0D0 

B(4,2)=0.0D0 

B(4,3)=0.0D0 

B(4,4)=1.0D0 

RETURN 

END 
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c 

c 

c 



c 



c 

c 



c 



c 



PROGRAM THESIS5.F0R (Time Delay-lst Order Approx in delta & y ie . in T1 & T2)) 
BIFURCATION ANALYSIS 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DOUBLE PRECISION K1 ,K2,L,NR,NV,IZ,MASS,NDELTA,NPSI ,NY, 

& NRDOT,NVDOT 

DIMENSION A(4,4) ,B(4,4) ,FV1(4) ,IV1 (4) , ZZZ(4,4) , ALFR(4) , 

& ALFI (4) ,BETA (4) ,WR(4) ,WI (4) 



OPEN (ll,FILE='BIFl.RESO 
OPEN (12,FILE='BIF2.RESO 
OPEN (13,FILE='BIF3.RESO 



Vehicle 


1 Parameters 


IZ 


=0.0 


L 


=528 


RHO 


= 1.94 


G 


=32.2 


XG 


=0.0 


MASS 


=0.0088 


U 


=3.0 



YRDOT 


= 0.00000 


YVDOT 


=-0.00912 


YR 


=+0.00456 


YV 


=-0.01434 


YPSI 


= 0.01400 


YY 


= 0.02000 



YDELTA= 0.00278 
NRDOT =-0.00115 
NVDOT = 0.00000 
NR =-0.00296 
NV =-0.00460 

NPSI = 0.01000 
NY = -0.00250 
NDELTA= -0.00139 



WRITE (*,1001) 
READ (*,*) 
WRITE (*,1002) 
READ (*,*) 
WRITE (*,1003) 
READ (*,*) 



WNMIN,WNMAX,IWN 

XDMIN,XDMAX,IXD 

ZETA 
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WRITE(*,1100) 

READ (♦,*) TLl 

TL1=TL1*U/L 

WRITE(*,1101) 

READ (*,*) TL2 

TL2=TL2*U/L 

DVR =(IZ-NRDOT)*(MASS-YVDOT)- 
& (MASS*XG-YRDOT) * (MASS*XG-NVDOT) 

AA11=( (IZ-NRDOT) *YPSI- (MASS*XG-YRDOT) *NPSI) /DVR 
AA12=( (IZ-NRDOT) *YV- (MASS*XG-YRDOT) *NV) /DVR 
AA21=((NVD0T-MASS*XG)*YPSI+(MASS-YVD0T)*NPSI)/DVR 
AA22= ( (MASS-YVDOT) *NV- (MASS*XG-NVDOT) *YV) /DVR 
AA13=( (IZ-NRDOT) *(YR-MASS)+ 

& (YRDOT-MASS*XG) * (NR-MASS*XG) ) /DVR 

AA23=( (NVDOT-MASS*XG) * (YR-MASS) + 

& (MASS-YVDOT) *(NR-MASS*XG)) /DVR 

AA 14= ( ( IZ-NRDOT) *YY+ (YRDOT-MASS*XG) *NY) /DVR 
AA24= ( (NVDOT-MASS*XG) *YY+ (MASS-YVDOT) *NY) /DVR 

BBl =( (IZ-NRDOT) *YDELTA- (MASS+XG-YRDOT) *NDELTA) /DVR 
BB2 = ( (MASS-YVDOT) *NDELTA- (MASS*XG-NVDOT) *YDELTA) /DVR 

AN0M=AA23 

BN0M=BB2 

CN0M=AA21 

EPS =l.D-5 

ILMAX=1500 

DO 1 1=1, IWN 

WRITE (*,2001) I, IWN 

WN=WNMIN+ ( I- 1 ) * ( WNMAX-WNMIN) / ( IWN- 1 ) 

Kl=- ( WN**2 . DO/BNOM) - (CNOM/BNOM) 

K2=- ( ANOM+2 . DO*ZETA*WN) /BNOM 
DO 2 J=1,IXD 

XD=XDMIN+ ( J- 1) ♦ (XDMAX-XDMIN) / (IXD-1) 

CALL LINEAR(K1,K2,XD,AA11,AA12,AA13,AA14, 

& AA2 1 , AA22 , AA23 , AA24 , BBl , BB2 , A , B , TLl , TL2) 

CALL RGG (4 , 4 , A , B , ALFR , ALFI , BETA , 0 , ZZZ , lERR) 

DO 11 IJE=1,4 

WR ( I JE) =ALFR ( I JE) /BETA ( IJE) 

WI ( IJE) =ALFI (IJE) /BETA ( IJE) 

11 CONTINUE 

CALL DSTABL(DEOS,WR,WI,FREQ) 
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IF (J.GT.l) GO TO 10 

DEOSOO=DEOS 

XDOO =XD 

LL=0 

GO TO 2 

10 DEOSNN=DEOS 
XDNN =XD 
PR=DEOSNN*DEOSOO 
IF (PR.GT.O.DO) GO TO 3 
LL=LL+1 

IF (LL.GT.3) STOP 1000 
IL=0 

XDO=XDOO 
XDN=XDNN 
DEOSO=DEOSOO 
DEOSN=DEOSNN 
6 XDL=XDO 
XDR=XDN 
DEOSL=DEOSO 
DEOSR=DEOSN 
XD=(XDL+XDR)/2.D0 

CALL LINEAR (K1 , K2 , XD , AAl 1 , AA12 , AA13 , AA14 , AA21 , AA22 , 
& AA23 , A A24 , BB 1 , BB2 , A , B , TL 1 , TL2 ) 

CALL RGG (4 , 4 , A , B , ALFR , ALFI , BETA , 0 , ZZZ , lERR) 

DO 12 IJE=1,4 

WR ( I JE ) = ALFR ( I JE) / BETA ( I JE ) 

WI ( I JE) =ALFI (I JE) /BETA ( I JE) 

12 CONTINUE 

CALL DSTABL(DEOS,WR,WI,FREQ) 

DEOSM=DEOS 

XDM=XD 

PRL=DEOSL*DEOSM 

PRR=DEOSR*DEOSM 

IF (PRL.GT.O.DO) GO TO 5 

XDO=XDL 

XDN=XDM 

DEOSO=DEOSL 

DEOSN=DEOSM 

IL=IL+1 

IF (IL.GT.ILMAX) STOP 3100 
DIF=DABS(XDL-XDM) 

IF (DIF.GT.EPS) GO TO 6 

XD=XDM 

GO TO 4 

5 IF (PRR.GT.O.DO) STOP 3200 
XDO=XDM 
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XDN=XDR 

DEOSO=DEOSM 

DEOSN=DEOSR 

IL=IL+1 

IF (IL.GT.ILMAX) STOP 3100 
DIF=DABS (XDM-XDR) 

IF (DIF.GT.EPS) GO TO 6 
XD=XDM 

4 LLL=10+LL 

WRITE (ILL,*) XD.WN 
3 XDOO=XDNN 

DEOSOO=DEOSNN 
2 CONTINUE 
1 CONTINUE 



1001 


FORMAT 


(’ 


ENTER 


MIN, 


MAX, AND INCREMENTS OF WN 


(dimensionless) ’ ) 


1002 


FORMAT 


(’ 


ENTER 


MIN, 


MAX, AND INCREMENTS OF XD 


(dimensionless) ’ ) 


1003 


FORMAT 


(’ 


ENTER 


DAMPING RATIO’) 




1100 


FORMAT 


(’ 


ENTER 


TIME 


LAG TLl (dimensional)’) 




1101 


FORMAT 


(’ 


ENTER 


TIME 


LAG TL2 (dimensional)’) 




2001 


FORMAT 


(215) 









END 



SUBROUTINE DSTABL (DEOS , WR , WI , OMEGA) 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DIMENSION WR(4),WI(4) 

DE0S=-1.0D+20 
DO 1 1=1,4 

IF (WR(I) .LT.DEOS) GO TO 1 
DEOS=WR(I) 

IJ=I 

1 CONTINUE 
OMEGA=WI(IJ) 

OMEGA=DABS (OMEGA) 

RETURN 

END 

SUBROUTINE LINEAR (K1 , K2 , XD , AAl 1 , AA12 , AA13 , AA14 , AA2 1 , AA22 , AA23 , AA24 
& ,BB1,BB2,A,B,TL1,TL2) 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DOUBLE PRECISION K1.K2 
DIMENSION A(4,4) ,B(4,4) 

A(1,1)=0.0D0 

A(1,2)=0.0D0 

A(1,3)=1.0D0 

A(1,4)=0.0D0 
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A(2,1)=AA11+BB1*K1-BB1*K1*TL2/XD-BB1*K1*TL1/XD 

A(2,2)=AA12-BB1*K1*TL2/XD-BB1*K1*TL1/XD 

A(2,3)=AA13+BB1*K2-BB1*K1*TL1+BB1*K1*TL1*TL2/XD 

A(2,4)=AA14+BB1*K1/XD 

A(3,1)=AA21+BB2*K1-BB2*K1*TL2/XD-BB2*K1*TL1/XD 

A(3,2)=AA22-BB2*K1*TL2/XD-BB2*K1*TL1/XD 

A(3,3)=AA23+BB2*K2-BB2*K1*TL1+BB2*K1*TL1*TL2/XD 

A(3,4)=AA24+BB2*K1/XD 

A(4,1)=1.0D0 

A(4,2)=1.0D0 

A(4,3)=0.0D0 

A(4,4)=0.0D0 

B(1,1)=1.0D0 

B(1,2)=O.ODO 

B(1,3)=O.ODO 

B(1,4)=0.0D0 

B(2,1)=0.0D0 

B (2 , 2)=1 . ODO- (BB1*K1*TL1*TL2/XD) 

B(2,3)=BB1*K2*TL1 

B(2,4)=0.0D0 

B(3,1)=0.0D0 

B(3,2)=-BB2*K1*TL1*TL2/XD 
B(3, 3)=1 . 0D0+(BB2*K2*TL1) 

B(3,4)=O.ODO 

B(4,1)=O.ODO 

B(4,2)=0.0D0 

B(4,3)=0.0D0 

B(4,4)=1.0D0 

RETURN 

END 
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PROGRAM HOPF.FOR 



C 
C 

C Hopf Bifurcation Analysis 

C 

C Third Order Expansions: First Order Approximation 

C 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DOUBLE PRECISION K1,K2,K3,L,NR,NV,NDELTA,NPSI,NY,IZ,MASS, 

& NRDOT , NVDOT , KIP , K2P , NWV , NVRR , NRVV , NYYY , NPPP 

DOUBLE PRECISION M11,M12,M13,M14,M21,M22,M23,M24, 

1 M31,M32,M33,M34,M41,M42,M43,M44, 

2 N11,N12,N13,N14,N21,N22,N23,N24, 

3 N31,N32,N33,N34,N41,N42,N43,N44, 

4 L21,L22,L23,L24,L31,L32,L33,L34, 

5 L41,L42,L43,L44 
C 

DIMENSION A(4,4) ,T(4,4) ,TINV(4,4) ,FV1(4) ,IV1(4) ,ZZZ(4,4) 
DIMENSION WR(4) ,WI(4) ,TSAVE(4,4) ,TLUD(4,4) , IVLUD(4) ,SVLUD(4) 
DIMENSION ASAVE(4,4) 



OPEN (11,FILE=^BIF1.RES0 
OPEN (15,FILE=^H0PF.RESO 
C 

C Vehicle Parameters: 



IZ 


=0.0 


L 


=528 


RHO 


=1.94 


G 


=32.2 


XG 


=0.0 


MASS 


=0.0088 


U 


=1.0 


YRDOT 


= 0.00000 


YVDOT 


=-0.00912 


YR 


=+0.00456 


YV 


=-0.01434 


YVW 


=-0.15391 


YVRR 


=-0.05476 


YRW 


= 0.04608 


YYYY 


= 0.46800 


YPPP 


= 0.00000 


YPSI 


= 0.01400 


YY 


= 0.02000 



YDELTA=+0. 00278 
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NRDOT =-0.00115 
NVDOT = 0.00000 
NR =-0.00296 
NV =-0 . 00460 
NVW =-0 . 00336 
NVRR = 0.00759 
NRVV =-0.05160 
NYYY = 0.00000 
NPPP = 0 . 00000 
NPSI = 0.01000 
NY =-0.00250 
NDELTA=-0. 00139 



WRITE (*,1003) 

READ (*,*) ZETA 

WRITE (*,1100) 

READ (*,*) TL 

TL=TL*U/L 
WRITE (*,1006) 

READ (*,*) DO 



C 

DVR =(IZ-NRDOT)*(MASS-YVDOT)- 
& (MASS*XG-YRDOT)*(MASS*XG-NVDOT) 

AA11=( (IZ-NRDOT) *YPSI-(MASS*XG-YRDOT) *NPSI)/DVR 
AA12= ( (IZ-NRDOT) *YV- (MASS*XG-YRDOT) *NV) /DVR 
AA21= ( (NVDOT-MASS*XG) *YPSI+ (MASS-YVDOT) *NPSI) /DVR 
AA22= ( (MASS-YVDOT) *NV- (MASS*XG-NVDOT) * YV) /DVR 
AA13= ( (IZ-NRDOT) * (YR-MASS) + 

& (YRDOT-MASS*XG) * (NR-MASS*XG) ) /DVR 

AA23= ( (NVDOT-MASS*XG) * (YR-MASS) + 

& (MASS-YVDOT) *(NR-MASS*XG)) /DVR 

AA14= ( (IZ-NRDOT) *YY+ (YRDOT-MASS*XG) *NY) /DVR 
AA24= ( (NVDOT-MASS*XG) *YY+ (MASS-YVDOT) *NY) /DVR 

BBl =( (IZ-NRDOT) *YDELTA- (MASS*XG-YRDOT) *NDELTA) /DVR 
BB2 = ( (MASS-YVDOT) *NDELTA- (MASS*XG-NVDOT) *YDELTA) /DVR 
C 

AN0M=AA23 

BN0M=BB2 

CN0M=AA21 

EPS =l.D-5 
ILMAX=1500 
C 

IWN=1000 
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DO 1 11=1, IWN 

WRITE (*,2001) II 

READdl,*) XD,WN 

Kl=- (WN**2 . DO/BNOM) - (CNOM/BNOM) 

K2=- ( ANOM+2 . DO*ZETA*WN) /BNOM 

K3=K2 



Start Hopf Bifurcation Analysis 

Evaluate Nonlinear Rudder Expansion Coefficients 



K2=0 . ODO 

K1P=K1-K1*TL*U/XD 

K2P=K2-K1*TL/XD 

DPPV=- ( 1 . DO/ (3 . D0*D0* *2) ) *3 . D0*K1P*K1P*K2P 

+ 0.5D0*K1*TL/XD + 3.D0*K1*(TL*U)**3/(3.D0*XD**3) 
DPVV=- ( 1 . DO/ (3 . D0*D0* *2) ) *3 . D0*K1P*K2P*K2P 
+ 3 . D0*U*TL*TL*TL*K1/ (3 . D0*XD**3) 

DPPR=- (1 . DO/ (3 . D0*D0**2) ) *3 . D0*K1P*K1P*K3 
DPRR=- ( 1 . DO/ (3 . D0*D0**2) ) *3 . D0*K1P*K3*K3 
DPPY=- (1 . DO/ (3 . D0*D0**2) ) *3 . D0*K1P*K1P*K1/XD 

- 3.D0*TL*TL*U*U*K1/(3.D0*XD**3) 

DPYY=- (1 . DO/ (3 . D0*D0**2) ) *3 . D0*K1P*K1*K1/ (XD**2) 

+ 3.D0*TL*U*K1/(3.D0*XD**3) 

DVVR=- (1 . DO/ (3 . D0*D0**2) ) *3 . D0*K2P*K2P*K3 
DVRR=- (1 . DO/ (3 . D0*D0**2) ) *3 . D0*K2P*K3*K3 
DVVY=- (1 . DO/ (3 . D0*D0**2) ) *3 . D0*K1*K2P*K2P/XD 

- 3.D0*K1*TL*TL/(3.D0*XD**3) 

DVYY=- ( 1 . DO/ (3 . D0*D0+*2) ) *3 . D0*K1*K1*K2P/ (XD**2) 

+ 3.D0*TL*K1/(3.D0*XD**3) 

DRRY=- (1 . DO/ (3 . D0*D0**2) ) *3 . D0*K1*K3*K3/XD 
DRYY=- (1 . DO/ (3 . D0*D0**2) ) *3 . D0*K1*K1*K3/ (XD**2) 

DPVR=- (1 . DO/ (3 . D0*D0**2) ) *6 . D0*K1P*K2P*K3 
DPVY=- (1 . DO/ (3 . D0*D0**2) ) *6 . D0*K1P*K1*K2P/XD 

- 6.D0*TL*TL*U*K1/(3.D0*XD**3) 

DPRY=- (1 . DO/ (3 . D0*D0**2) ) *6 . D0*K1P*K1*K3/XD 
DVRY=- ( 1 . DO/ (3 . D0*D0**2) ) *6 . D0*K1*K2P*K3/XD 
DPPP=- (1 . DO/ (3 . D0*D0**2) ) *1 . D0*K1P*K1P*K1P 

+ K1*TL*U/(6.D0*XD) + (K1*(TL*U)**3)/(3.D0*XD**3) 
DVVV=- (1 .DO/ (3 . D0*D0**2) ) *1 .D0*K2P*K2P*K2P 
+ K1*(TL**3)/(3.D0*XD**3) 

DRRR=- (1 . DO/ (3 . D0*D0**2) ) *1 . D0*K3*K3*K3 

DYYY=- (1 . DO/ (3 . D0*D0**2) ) * (Kl/XD) **3-Kl/ (3 . D0*XD**3) 

- K1/(3.D0*XD**3) 



Evaluate Transformation Matrix of Eigenvectors 



CALL LINEAR (K1 , K3 , XD , AAl 1 , AA12 , AA13 , AA 14 , AA21 , AA22 , AA23 , 
& AA24,BB1,BB2,A,TL) 

DO 11 1=1,4 
DO 12 J=l,4 

ASAVE(I,J)=A(I, J) 

12 CONTINUE 
11 CONTINUE 

CALL RG (4 . 4 , A . WR . WI . 1 , ZZZ . I VI . FVl . lERR) 

CALL DSOMEG(IEV,WR,WI, OMEGA, CHECK) 

OMEGAO=OMEGA 
DO 50 1=1,4 

T(I,1)= ZZZ(I,IEV) 

T(I,2)=-ZZZ(I,IEV+1) 

50 CONTINUE 

IF (lEV.EQ.l) GO TO 13 

IF (IEV.EQ.2) GO TO 14 

IF (IEV.EQ.3) GO TO 15 

STOP 3004 

14 DO 60 1=1,4 

T(I,3)=ZZZ(I,1) 

T(I,4)=ZZZ(I,4) 

60 CONTINUE 
GO TO 17 

15 DO 70 1=1,4 

T(I,3)=ZZZ(I,1) 

T(I,4)=ZZZ(I,2) 

70 CONTINUE 
GO TO 17 

13 DO 16 1=1,4 

T(I,3)=ZZZ(I,3) 

T(I,4)=ZZZ(I,4) 

16 CONTINUE 

17 CONTINUE 
C 

C Normalization of the Critical Eigenvector 

C 

CALL NORMAL (T) 

C 

C Definition of Mij 

C 

M11=T(1,1) 

M21=T(2,1) 

M31=T(3,1) 

M41=T(4,1) 

M12=T(1,2) 

M22=T(2,2) 
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c 

c 

c 



M32=T(3,2) 

M42=T(4,2) 

M13=T(1,3) 

M23=T(2,3) 

M33=T(3,3) 

M43=T(4,3) 

M14=T(1,4) 

M24=T(2,4) 

M34=T(3,4) 

M44=T(4,4) 



Definition of Lij 



L21= DPPV+M11+M11+M21 
& + DPRR*M11*M31*M31 
& + DVVR*M21*M21*M31 
& + DVYY+M21+M41+M41 
& + DPVR+M11+M21+M31 
& + DVRY*M21*M41*M31 
& + DRRR*M31*M31*M31 



+ DPVV*M11*M21*M21 
+ DPPY*M11*M11*M41 
+ DVRR+M21+M31+M31 
+ DRRY+M31+M31+M41 
+ DPVY*M11*M21*M41 
+ DPPP*M11*M11*M11 
+ DYYY*M41*M41*M41 



+ DPPR*M11*M11*M31 
+ DPYY*M11*M41*M41 
+ DVVY*M21*M21*M41 
+ DRYY+M31+M41+M41 
+ DPRY+M11+M41+M31 
+ DVVV*M21*M21*M21 



PPV = M11*M11*M22 + 2.0*M11*M12*M21 
PVV = M12*M21*M21 + 2. 0*M11*M21*M22 
PPR = M11+M11+M32 + 2.0+M11+M12+M31 
PRR = M12*M31*M31 + 2 . 0*M11*M31*M32 
PPY = M11*M11*M42 + 2.0*M11*M12*M41 
PYY = M41*M41*M12 + 2. 0*M11*M41*M42 
WR = M21+M21+M32 + 2 . 0*M31*M21*M22 
VRR = M22+M31+M31 + 2 . 0+M31+M32+M21 
VVY = M21*M21*M42 + 2 . 0+M41+M21+M22 
VYY = M22*M41*M41 + 2 . 0*M41*M42*M21 
RRY = M31*M31*M42 + 2 . 0*M41*M31*M32 
RYY = M32*M41*M41 + 2 . 0*M41*M42*M31 
PVR = M11*M21*M32+M31*(M11*M22+M12*M21) 

PVY = M11*M21*M42+M41*(M11*M22+M12*M21) 

PRY = M11*M41*M32+M31*(M11*M42+M12*M41) 

VRY = M21*M41*M32+M31*(M21*M42+M22*M41) 

PPP = 3.0+M11+M11+M12 
VVV = 3.0*M21*M21*M22 
RRR = 3.0*M31*M31*M32 
YYY = 3.0*M41*M41*M42 

L22=DPPV*PPV+DPVV*PW+DPPR*PPR+DPRR*PRR+DPPY*PPY+DPYY*PYY 

+DWR*VVR+DVRR*VRR+DVVY*VVY+DVYY*VYY+DRRY*RRY+DRYY*RYY 

+DPVR*PVR+DPVY*PVY+DPRY*PRY+DVRY*VRY+DPPP*PPP+DVW*VVV 

+DRRR+RRR+DYYY+YYY 



PPV = M12+M12+M21 + 2.0*M11*M12*M22 
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c 



c 



c 



c 



PVV = M11*M22*M22 + 2. 0*M12*M21*M22 
PPR = M12*M12*M31 + 2. 0*M11*M12*M32 
PRR = M11*M32*M32 + 2. 0*M12*M31*M32 
PPY = M12*M12*M41 + 2. 0*M11*M12*M42 
PYY = M11*M42*M42 + 2 . 0*M12*M41*M42 
WR = M22*M22*M31 + 2 . 0*M21*M22*M32 
VRR = M21*M32*M32 + 2. 0*M22*M31*M32 
VVY = M22*M22*M41 + 2 . 0*M21*M22*M42 
VYY = M21*M42*M42 + 2 . 0*M22*M41*M42 
RRY = M32*M32*M41 + 2. 0*M31*M32*M42 
RYY = M31*M42*M42 + 2 . 0*M32*M41*M42 
PVR = M12*M22*M31 + M32* (M11*M22+M12*M21) 

PVY = M12*M22*M41 + M42* (M11*M22+M12*M21) 

PRY = M12*M42*M31 + M32* (M11*M42+M12*M41) 

VRY = M22*M42*M31 + M32* (M21*M42+M22*M41) 

PPP = 3.0*M11*M12*M12 
WV = 3.0*M21*M22*M22 
RRR = 3.0*M31*M32*M32 
YYY = 3.0*M41*M42*M42 

L23=DPPV*PPV+DPVV*PW+DPPR*PPR+DPRR*PRR+DPPY*PPY+DPYY*PYY 
& +DWR*VVR+DVRR*VRR+DVVY*WY+DVYY*VYY+DRRY*RRY+DRYY*RYY 

& +DPVR*PVR+DPVY*PVY+DPRY*PRY+DVRY*VRY+DPPP*PPP+DVW*VVV 

& +DRRR*RRR+DYYY*YYY 



& 

& 

& 

& 

& 

& 



L24= DPPV*M12*M12*M22 
+ DPRR*M12*M32*M32 
+ DVVR*M22*M22*M32 
+ DVYY*M22*M42*M42 
+ DPVR*M12*M22*M32 
+ DVRY*M22*M42*M32 
+ DRRR*M32*M32*M32 



+ DPVV*M12*M22*M22 
+ DPPY*M12*M12*M42 
+ DVRR*M22*M32*M32 
+ DRRY*M32*M32*M42 
+ DPVY*M12*M22*M42 
+ DPPP*M12*M12*M12 
+ DYYY*M42*M42*M42 



+ DPPR*M12*M12*M32 
+ DPYY*M12*M42*M42 
+ DVVY*M22*M22*M42 
+ DRYY*M32*M42*M42 
+ DPRY*M12*M42*M32 
+ DVVV*M22*M22*M22 



L31=L21 

L32=L22 

L33=L23 

L34=L24 



L21=L21*BB1*U*U 

L22=L22*BB1*U*U 

L23=L23*BB1*U*U 

L24=L24*BB1*U*U 

L31=L31*BB2*U*U 

L32=L32*BB2*U*U 

L33=L33*BB2*U*U 

L34=L34*BB2*U*U 

L41=-0 . 1* (M21+U+M1 1/3 . 0) 

L42=-M11* (M12+M21+0 . 5*M1 l*M22+0 . 5*U*M12*M11) 
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L43=-M12* (Ml 1 *M22+0 . 5*M12*M21+0 . 5*U*M11*M12) 

L44=-0. 5*M12*M12* (M22+U*M12/3 . 0) 

C 

Cl 1= ( 1/6 . 0) *YVW* (M21*M21*M21)+0 . 5*YVRR* (M31*M31 ♦M21) + 

& 0.5*YRVV*M21*M21*M31+(1/6.0)*YPPP*M11*M11*M11+(1/6.0)*YYYY*M41*M41*M41 

C12= ( 1/6 . 0) ♦YVW* (3*M21*M21*M22)+0 . 5*YVRR* (M31*M31* 

& M22+2*M31*M32*M21)+0.5*YRVV*(M21*M21*M32+2*M31*M21*M22)+(1/6.0)*YPPP*3*M11*M11*M12 

& +(1/6.0)*YYYY*3*M41*M41*M42 



C13= ( 1/6 . 0) *YVW* (3*M21*M22*M22)+0 . 5*YVRR* (M21*M32* 

& M32+2*M31*M32*M22)+0 . 5*YRW* (M22*M22*M31+2*M32*M21*M22)+(l/6 . 0) *YPPP*3*M11*M12*M12 

& +(1/6.0)*YYYY*3*M41*M42*M42 

C14= ( 1/6 . 0) *YVW* (M22*M22*M22)+0 . 5*YVRR* (M32*M32*M22) + 

& 0.5*YRVV*M22*M22*M32+(1/6.0)*YPPP*M12*M12*M12+(1/6.0)*YYYY*M42*M42*M42 

C21= ( 1/6 . 0) *NVW* (M21*M21*M21)+0 . 5* (NVRR*M31*M31*M21) + 

& 0.5*NRVV*M21*M21*M31+(1/6.0)*NPPP*M11*M11*M11+(1/6.0)*NYYY*M41*M41*M41 

C22= (1/6 . 0) ♦NVW* (3*M21*M21*M22)+0 . 5*NVRR* (M31*M31* 

& M22+2*M31*M32*M21)+0.5*NRW*(M21*M21*M32+2*M31*M21*M22)+(1/6.0)*NPPP*3*M11*M11*M12 

& +(1/6.0)*NYYY*3*M41*M41*M42 

C23= ( 1/6 . 0) ♦NVW* (3*M21*M22*M22)+0 . 5*NVRR* (M21*M32* 

& M32+2*M31*M32*M22)+0.5*NRW*(M22*M22*M31+2*M32*M21*M22)+(1/6.0)*NPPP*3*M11+M12*M12 

& +(1/6.0)*NYYY*3*M41*M42*M42 

C24= (1/6.0) *NVW* (M22*M22*M22) +0 . 5*NVRR*M32*M32*M22+ 

& 0.5*NRW*M22*M22*M32+(1/6.0)*NPPP*M12*M12*M12+(1/6.0)*NYYY*M42*M42*M42 



D11=(C11*(IZ-NRD0T)+C21*(YRD0T-M*XG))/DVR 
D12=(C12* (IZ-NRD0T)+C22* (YRD0T-M*XG) ) /DVR 
D13= (C13* ( IZ-NRDOT) +C23* (YRD0T-M*XG) ) /DVR 
D14= (C14* (IZ-NRDOT) +C24* (YRDOT-M*XG) ) /DVR 

D21= (Cl 1* (NVD0T-M*XG)+C21* (M-YVDOT) ) /DVR 
D22=(C12*(NVD0T-M*XG)+C22*(M-YVD0T))/DVR 
D23= (C13* (NVDOT-M*XG) +C23* (M-YVDOT) ) /DVR 
D24= (C14* (NVD0T-M*XG)+C24* (M-YVDOT) ) /DVR 
C 

C Invert Transformation Matrix 

C 

DO 20 1=1,4 
DO 30 J=l,4 
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TINV(I,J)=0.0 

TSAVE(I,J)=T(I,J) 

30 CONTINUE 

20 CONTINUE 

CALL DLUD(4,4,TSAVE,4,TLUD,IVLUD) 

DO 40 1=1,4 

IF (IVLUD(I) .EQ.O) STOP 3003 
40 CONTINUE 

CALL DILU(4,4,TLUD,IVLUD,SVLUD) 

DO 80 1=1,4 
DO 90 J=l,4 

TINV(I,J)=TLUD(I, J) 

90 CONTINUE 

80 CONTINUE 
C 

C Check Inv(T)*A*T 

C 

IMULT=2 

IF (IMULT.EQ.l) CALL MULT(TINV, ASAVE,T) 

IF (IMULT.EQ.O) STOP 
C 

C Definition of Nij 

C 

N11=TINV(1,1) 

N21=TINV(2,1) 

N31=TINV(3,1) 

N41=TINV(4,1) 

N12=TINV(1,2) 

N22=TINV(2,2) 

N32=TINV(3,2) 

N42=TINV(4,2) 

N13=TINV(1,3) 

N23=TINV(2,3) 

N33=TINV(3,3) 

N43=TINV(4,3) 

N14=TINV(1,4) 

N24=TINV(2,4) 

N34=TINV(3,4) 

N44=TINV(4,4) 

C 

R11=N12*L21+N13*L31+N14*L41+N12*D11+N13*D21 

R12=N12*L22+N13*L32+N14*L42+N12*D12+N13*D22 

R13=N12*L23+N13*L33+N14*L43+N12*D13+N13*D23 

R14=N12*L24+N13*L34+N14*L44+N12*D14+N13*D24 

R21=N22*L21+N23*L31+N24*L41+N22*D11+N23*D21 

R22=N22*L22+N23*L32+N24*L42+N22*D12+N23*D22 

R23=N22*L23+N23*L33+N24*L43+N22*D13+N23*D23 

R24=N22*L24+N23*L34+N24*L44+N22*D14+N23*D24 
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c 

C Evaluate Alpha ^ and Omega ^ 

C 

DINC=0.001 
XDR =XD+DINC 
XDL =XD-DINC 
XD =XDR 

CALL LINEAR (K1 , K3 , XD , AAl 1 , AA12 , AA13 , AAl 4 , AA21 , AA22 , AA23 , 
& AA24,BB1,BB2,A,TL) 

CALL RG(4,4,A,VR,WI,0,ZZZ,IV1,FV1,IERR) 

CALL DSTABL(DEOS,WR,WI,FREQ) 

ALPHR=DEOS 

OMEGR=FREQ 

XD=XDL 

C 

CALL LINEAR (K1 , K3 , XD , AAl 1 , AA12 , AA13 , AA14 , AA21 , AA22 , AA23 , 
& AA24,BB1,BB2,A,TL) 

CALL RG(4,4,A,WR,WI,0,ZZZ,IV1,FV1,IERR) 

CALL DSTABL(DEOS,WR,WI,FREQ) 

ALPHL=DEOS 

OMEGL=FREQ 

C 

DALPHA= ( ALPHR-ALPHL) / (XDR-XDL) 

DOMEGA= (OMEGR-OMEGL) / (XDR-XDL) 

C 

C Evaluation of Hopf Bifurcation Coefficients 

C 

COEF 1 =3 . 0*R1 1+R13+R22+3 . 0*R24 
C0EF2=3 . 0*R21+R23-R12-3 . 0*R14 
AMU2 =-C0EFl/(8.0*DALPHA) 

BETA2=0.25*C0EF1 

TAU2 =- (C0EF2-D0MEGA*C0EF1/DALPHA) / (8 . 0*OMEGAO) 

PER =2 . 0*3 . 1415927/OMEGAO 
PER =PER*U/L 
C 

WRITE (15,2002) XD,WN,C0EF1 ,DALPHA,PER 



1 CONTINUE 
C 

1001 FORMAT C ENTER MIN, MAX, AND INCREMENTS OF WN (dimensionless) O 

1002 FORMAT (' ENTER MIN, MAX, AND INCREMENTS OF XD (dimensionless) O 

1003 FORMAT (’ ENTER DAMPING RATIO’) 
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1006 FORMAT (’ ENTER DSAT ’) 

1100 FORMAT (' ENTER TIME LAG TL (dimensional)’) 

2001 FORMAT (215) 

2002 FORMAT (5E15.5) 

END 

SUBROUTINE DSTABL(DEOS,WR,WI, OMEGA) 

C 

C EVALUATES THE EIGENVALUE WITH THE MAXIMUM REAL PART 
C 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

DIMENSION WR(4),WI(4) 

DE0S=-1.0D+20 
DO 1 1=1,4 

IF (WR(I) .LT.DEOS) GO TO 1 
DEOS=WR(I) 

U=I 

1 CONTINUE 
OMEGA=WI(IJ) 

OMEGA=DABS (OMEGA) 

RETURN 

END 

C 

SUBROUTINE LINEAR (K 1 , K2 , XD , AAl 1 , AA12 , AA13 , AA14 , AA2 1 , 

& AA22,AA23,AA24,BB1,BB2,A,TL) 

C 

C FORMS THE LINEARIZED MATRIX A (time delay 1st order approximation 
C 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DOUBLE PRECISION K1,K2 
DIMENSION A(4,4) 

A(1,1)=O.ODO 

A(1,2)=0.0D0 

A(1,3)=1.0D0 

A(1,4)=0.0D0 

A(2,1)=AA11+BB1*K1-BB1*K1*TL/XD 

A (2 , 2) =AA12-BB1*K1*TL/XD 

A(2,3)=AA13+BB1*K2 

A(2,4)=AA14+BB1+K1/XD 

A(3,1)=AA21+BB2*K1-BB2*K1*TL/XD 

A (3 , 2) =AA22-BB2*K1*TL/XD 

A(3,3)=AA23+BB2*K2 

A(3,4)=AA24+BB2*K1/XD 

A(4,1)=1.0D0 

A(4,2)=1.0D0 



79 



A(4,3)=0.0D0 

A(4,4)=0.0D0 

RETURN 

END 



C 



SUBROUTINE DSOMEG ( UK , WR , WI , OMEGA , CHECK) 
IMPLICIT DOUBLE PRECISION (A-H.O-Z) 
DIMENSION WR(4),WI(4) 

CHECK=-1.0D+25 
DO 1 1=1,4 

IF (WR(I) .LT.CHECK) GO TO 1 
CHECK=WR(I) 

U=I 

1 CONTINUE 

OMEGA=DABS(WI(IJ)) 

IF (WI(IJ) .GT.O.DO) IJK=IJ 
IF (WI(IJ) .LT.O.DO) IJK=IJ-1 
RETURN 
END 



SUBROUTINE NORMAL (T) 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

DIMENSION T(4,4) ,TN0R(4,4) 
CFAC=T(1,1)**2+T(1,2)**2 
IF (DABS(CFAC) .LE. (l.D-10)) STOP 4001 
TN0R(1,1)=1.D0 

TN0R(2,1)=(T(1,1)*T(2,1)+T(2,2)*T(1,2))/CFAC 

TN0R(3,1)=(T(1,1)*T(3,1)+T(3,2)*T(1,2))/CFAC 

TN0R(4,1)=(T(1,1)*T(4,1)+T(4,2)*T(1,2))/CFAC 

TNOR(1,2)=O.DO 

TN0R(2,2)=(T(2,2)*T(1,1)-T(2,1)*T(1,2))/CFAC 
TN0R(3,2)=(T(3,2)*T(1,1)-T(3,1)*T(1,2))/CFAC 
TN0R(4,2)=(T(4,2)*T(1,1)-T(4,1)*T(1,2))/CFAC 
DO 1 1=1,4 
DO 2 J=l,2 

T(I, J)=TNOR(I, J) 

2 CONTINUE 
1 CONTINUE 
RETURN 
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END 



SUBROUTINE MULT(TINV, A ,T) 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DIMENSION TINV(4,4) ,A(4,4) ,T(4,4) ,A1(4,4) ,A2(4,4) 
DO 1 1=1,4 
DO 2 J=l,4 
Aid, J)=O.DO 
A2(I, J)=O.DO 

2 CONTINUE 
1 CONTINUE 

DO 3 1=1,4 
DO 4 J=l,4 
DO 5 K=l,4 

Aid, J)=A(I,K)*T(K,J)+A1(I, J) 

5 CONTINUE 
4 CONTINUE 

3 CONTINUE 
DO 6 1=1,4 

DO 7 J=l,4 
DO 8 K=l,4 

A2(I , J)=TINV(I , K) *A1 (K, J)+A2(I , J) 

8 CONTINUE 

7 CONTINUE 

6 CONTINUE 

DO 11 1=1,4 

C WRITE (*,101) (Ad, J) , J=l,4) 

11 CONTINUE 

DO 12 1=1,4 

C WRITE (*,101) (T(I, J) , J=l,4) 

12 CONTINUE 

DO 10 1=1,4 

WRITE (*,101) (A2(I, J) , J=l,4) 

10 CONTINUE 

C WRITE (*,101) A2(l,l) 

RETURN 

101 FORMAT (4E15.5) 

END 
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